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BLOCK  19:  Abstract 


The  second  and  third  papers  are  on  the  subject  of  sources  embedded  in  anisotropic  media.  In  the 
second  paper  asymptotic  far  field  expressions  of  the  exact  spectral  elastodynamic  Green’s  tensor  are  found; 
in  the  third  paper  this  tensor  is  evaluated  numerically  by  the  discrete  wavenumber  method.  These  results 
are  used  to  examine  the  radiation  pattern  of  seismic  waves  from  different  types  of  sources.  This  is  a  crucial 
problem  in  monitoring  underground  nuclear  tests  because  the  radiation  pattern  is  a  potential  discriminant 
between  explosions  and  other  types  of  source.  Specifically,  it  is  demonstrated  that  for  specific  values  of 
the  parameters  of  the  anisotropic  medium  an  explosion  may  show  a  radiation  pattern  similar  to  that  of 
an  earthquake  in  an  isotropic  medium.  This  is  the  first  time  the  problem  of  a  source  embedded  in  an 
anisotropic  medium  has  been  solved;  the  effects  discussed  here  are  totally  separate  from  the  effects  of  plane 
wave  propagation  in  anisotropic  media  discussed  in  the  previous  literature. 

The  final  paper  treats  the  problem  of  attenuation  of  shear  wave  energy  in  the  crust.  Models  of  crustal 
attenuation  (Q)  as  a  function  of  depth  and  frequency  are  found  for  two  non-tectonic  regions,  the  north¬ 
eastern  U.S.  and  Fennoscandia,  with  crystalline  basement  rock  extending  to  the  surface.  The  method  used 
is  matching  complete  synthetic  seismograms  with  observed  seismograms.  Our  study  differs  from  other 
studies,  however,  in  that  we  are  interested  in  finding  a  general  model  of  crustal  attenuation  in  terms  of  the 
operative  processes.  Such  a  model  would  in  principle  be  applicable  to  any  geologically  similar  region.  We 
find  that  attenuation  is  high  ( Q  <  100  at  1  Hz)  in  the  upper  crust  and  is  probably  controlled  by  fluid  flow  in 
fractures  (frequency  dependence  of  Q  /°  5).  In  the  middle  crust  (depth  2-10  km)  attenuation  is  moderate 
(100  <  Q  <  500  at  1  Hz)  and  may  be  due  to  scattering  (frequency  dependence  of  Q  fl).  Attenuation 
is  low  and  frequency-independent  in  the  lower  crust,  indicating  anelastic  attenuation  in  an  environment 
where  microfractures  have  been  annealed. 
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PERTURBATION  APPROXIMATION  OF  3-D  SEISMIC  SCATTERING 


M.  P range  and  M.N.  Toksoz 
Earth  Resources  Laboratory 

Department  of  Earth,  Atmospheric,  and  Planetary  Sciences 
Massachusetts  Institute  of  Technology 

ABSTRACT 

A  method  is  presented  for  computing  three-dimensional  seismic  wave  scattering  from  a  rough  in¬ 
terface.  The  matrix  method  used  is  appropriate  for  direct  implementation  in  existing  propagator 
matrix-based  seismogram  synthesis  programs.  It  is  derived  using  a  perturbation  approach  which 
requires  interface  height  perturbations  to  be  small  relative  to  the  wavelengths  of  scattered  waves, 
and  interface  slope  perturbations  to  be  much  less  than  unity.  These  validity  conditions  are  based 
on  an  order-of-error  analysis  of  the  truncation  of  the  perturbation  series.  These  conditions  are 
numerically  investigated  by  comparison  of  frequency-wavenumber  domain  and  time  domain  per¬ 
turbation  results  with  those  generated  by  a  second-order  finite  difference  method  for  several  rough 
interface  models  with  Gaussian  auto-correlation  functions.  In  the  u-k  domain  comparisons,  the 
perturbation  method  is  accurate  for  RMS  interface  height  deviations  of  less  than  about  10  percent 
of  the  smallest  wavelength  in  the  scattered  field.  This  result  is  independent  of  RMS  interface  slope 
in  the  tested  range  of  0.037  to  0.99.  Comparisons  of  seismograms  generated  by  the  two  methods 
show  that  error  does  increase  with  increasing  RMS  slope,  but  at  half  the  rate  of  error  growth  with 
increasing  height.  Time  domain  error  is  acceptable  for  RMS  height  deviations  of  less  that  about  20 
percent  and  RMS  slopes  of  less  than  about  0.25.  A  three-dimensional  scattering  kernel  is  defined 
which  facilitates  analysis  of  two-  and  three-dimensional  scattered  field  results. 

INTRODUCTION 

The  presence  of  a  rough  interface  can  strongly  affect  seismic  waves  reflected  from  and  transmitted 
through  that  interface,  even  when  the  scale  of  roughness  is  much  less  than  a  wavelength.  These 
effects  include  changes  in  the  amplitude,  scattering  angle,  frequency  content,  and  wave-type  conver¬ 
sion  of  the  scattered  wave.  Available  exact  solutions  take  the  form  of  integral  equations  (DeSanto 
and  Brown,  198G)  or  finite  difTcrence/finite  element  formulations  (Levander  and  Hill,  1985)  that 
are  prohibitively  expensive  to  solve  in  three  dimensions.  In  this  paper,  we  present  a  perturba- 
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tion  approach  to  the  solution  of  the  three-dimensional  elastic  wave  equation  which  satisfies  welded 
boundary  conditions  at  a  rough  interface.  This  solution  is  an  extension  of  two-dimensional  surface 
wave  scattering  formulations  (Kennett,  1972;  Gilbert  and  Knopoff,  1960).  The  method  requires  the 
height  and  slope  of  interface  irregularities  to  be  small  with  respect  to  the  wavelengths  of  the  elas¬ 
tic  waves  present.  The  accuracy  of  the  perturbation  method  is  then  explored  for  two-dimensional 
models  by  comparing  solutions  for  a  series  of  rough  interface  models  with  those  generated  by  a 
second-order  finite  difference  method.  The  scattering  results  are  first  compared  in  the  frequency- 
wavenumber  domain  in  the  form  of  scattering  coefficients.  The  scattered  wave  fields  computed  with 
the  perturbation  and  finite  difference  methods  are  separated  into  up-  and  down-going  P  and  SV 
waves,  and  these  scattering  coefficients  are  individually  compared.  Comparisons  are  also  made  in 
the  time-space  domain  by  comparing  seismograms.  Finally,  we  discuss  the  features  of  the  three- 
dimensional  scattered  field. 

THREE-DIMENSIONAL  SCATTERING  FORMULATION 

A  fine-scale  blow-up  of  a  three-dimensional  rough  interface  is  shown  in  Figure  1.  The  irregular 
interface  is  described  by  z  —  h(x,y),  and  has  a  downward  normal  n(x,y).  It  separates  an  upper 
medium,  described  by  compressional  and  shear  wave  velocities  a i  and  /3y  and  density  py,  from  a 
lower  medium,  described  by  q2,  and  pi.  The  essence  of  the  formulation  is  to  project  displace¬ 
ments  and  stresses  on  the  two  sides  of  the  rough  interface  onto  a  planar  surface  whose  depth  equals 
the  mean  depth  of  the  rough  interface.  These  projected  fields  are  then  expanded  in  a  perturba¬ 
tion  series  in  h  about  a  background  field  consisting  of  the  known  planar  interface  solution.  This 
procedure  results  in  a  formulation  in  which  the  rough  interface  scattering  problem  is  replaced  by 
a  planar  interface  scattering  problem  with  sources  along  the  planar  interface  generating  the  rough 
interface  component  of  the  scattered  field. 

The  general  form  of  the  elastic  wave  equation  which  is  valid  for  small  displacements  in  the 
absence  of  body  forces  is 

-  pcv'2u,  =  Tjij.  i,j  e  (1) 

where  it,  is  the  ith  component  of  the  displacement  vector,  rJt  is  the  i,j  component  of  the  Cauchy 
stress  tensor,  the  comma  denotes  differentiation  of  tjx  in  the  j-th  coordinate  direction,  x  is  the 
temporal  frequency,  and  the  Einstein  summation  convention  applies.  Throughout  this  paper,  the 
Fourier  transform  in  x,  y,  and  t  uses  an  implied  phase  factor  of  exp (ikrx  +  ik^y  -  i~t).  For  an 


isotropic  solid,  stresses  are  related  to  displacements  by  the  constitutive  relation 


Ttj  =  \ukikf>i:  +  (2) 

where  6tJ  is  the  Kroniker  symbol,  and  A  and  \i  are  the  Lame  parameters.  Equations  (1)  and  (2) 
can  be  expressed  as  a  6  x  6  matrix  wave  equation  in  the  form 

'^  =  Ar ,  (3) 

dz  — 

where  r  is  the  displacement-stress  vector  defined  by  r  =  [ux,  uy,  uz,  txz ,  ryz,  ^2]r  and  ^  is 
defined  by 
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with  (  =  4/r(A  -f  /r)/(A  -f  2/z). 

Welded  boundary  conditions  at  the  rough  interface  require  continuity  of  displacement  and  trac¬ 
tion  at  each  point  on  the  interface.  These  tractions  are  measured  with  respect  to  the  local  tangent 
plane  at  each  point  on  the  interface,  and  are  given  by  T:  =  <7]knk,  where  the  downward  unit  normal 
nk  is  defined  by 


'1  +  h%  +  h\ 


A  new  displacement-stress  vector  f  is  defined  using  these  tangent  plene  tractions  so  that  r  is 
continuous  at  the  rough  interface,  r  is  defined  as 
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(6) 

Writing  (6)  in  a  form  which  explicitly  expresses  the  h  dependence,  this  transformation  takes  the 
form 


r{x,y\h)  =  (/  +  /i,x£x  +  h,yQy)r(x,y,h)  (7) 

where  r  and  r  are  evaluated  along  the  rough  interface,  £  is  the  identity  matrix,  and  Q  x  and  Q_y 
are 
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The  rough  interface  boundary  conditions  may  then  be  expressed  in  the  form 


L  (1)(x.y;/i)  =  f  (2)(x,y;/i)  (10) 

where  superscripts  indicate  the  respective  media. 

In  order  to  relate  the  scattered  field  (the  rough  interface  solution)  to  the  background  field  (the 
mean  planar  interface  solution),  the  scattered  field  at  the  mean  planar  interface  r  (0)  is  extrapolated 
to  the  rough  interface  by  the  power  series  expansion 

/i2 

r(h)  =  r( 0)  +  hr^{0)  +  -Tj-r.^O)  +  • . .  (11) 

Makii  g  use  of  the  wave  equation  (3),  (11)  is  reduced  to  the  form 

r(h)  =  {I  +  h±  +  ^A2  +  ...)L(0)  .  (12) 

So  far  the  formulation  is  exact  as  long  as  the  series  in  (12)  converges.  It  is  easy  to  demonstrate 
convergence  for  the  case  of  a  planar  interface.  In  this  case  h  is  constant  and  (12)  can  be  easily 
Fourier  transformed  to  the  (kx,ky)  domain  to  form  the  power  series  expansion  of  the  exponential 
function,  converging  to 


r(kx,ky-h )  =  eh=r(kx.ky-,0)  .  (13) 

This  is  the  standard  form  for  the  propagator  matrix  (Aki  and  Richards,  1980.  p.  275),  which  is  an 
exact  extrapolation  operator  which  forms  the  basis  of  the  propagator  matrix  method  for  formulating 
wave  equation  solutions  in  plane  layered  media  (Kennett,  1983b 

After  the  displacement-stress  vectors  along  the  rough  interface  have  been  extrapolated  to  the 
mean  planar  interface,  they  are  expanded  in  a  perturbation  series  about  ro(0),  the  d’splacement- 
stress  vector  at  the  interface  of  the  background  planar  interface  model: 

rW(0)  =  Zlo(0)  +  hr  [j)(0)  +  h2r[j)(  0)  +  . . .  ,  (14) 

where  superscripts  indicate  the  respective  media  and  subscripts  indicate  the  approximation  order  of 
the  field.  Since  ro(0)  is  the  displacement-stress  vector  for  the  background  planar  interface  model, 
it  is  continuous  across  the  planar  interface  and  needs  no  superscript.  The  higher  order  terms  reflect 
the  influence  of  the  rough  interface.  Combining  (7),  (12).  and  (14),  each  side  of  the  boundary 
condition  expressed  in  (10)  can  be  written  as  (omitting  the  superscripts) 


(15) 


r(h)  =  (/  +  /i,i£r  +  h,yQ_  y) 

(I  +  +  '2[42  +  •  •  ■) 

■  (lo  4-  hr  I  +  h2r2  +  . . .) 

~  (Z  +  ^,i£x  +  +  /l-‘  (16) 

+  0(h2)  +  0(hh'X)  +  0{hh,y) 

where  O(-)  denotes  order  of  accuracy. 

Applying  the  boundary  condition  (10)  to  (16)  results  in 


Ml(i2) -Li1')  ~  «(4(1) -4(2,)r0  +  M£y) -£x2))z:o  (17) 

+  h<y(gW-gW)r0, 

which  is  accurate  to  second-order  in  h  and  its  derivatives.  The  right-hand  side  of  (17)  is  clearly  zero 
in  a  planar  interface  model.  The  presence  of  interface  roughness  results  in  discontinuities  in  the 
displacement-stress  vector.  Such  discontinuities  represent  sources(Aki  and  Richards,  1980.  p.  38). 
Thus,  to  first  order  in  h,  the  effects  of  a  rough  interface  can  be  duplicated  by  adding  sources  along 
the  mean  planar  interface.  These  sources  will  be  designated  by  s  =  /t(r  [2)  -  rj1').  This  use  of 
sources  to  represent  material  deviations  from  a  background  model  is  shared  with  standard  Born 
theoretical  developments( Wu  and  Aki,  1985).  The  mapping  of  heterogeneity  into  source  terms  is 
also  used  in  exact  formulations  based  on  Huygen’s  principle  (Paul  and  Campillo,  1988). 

Fourier  transforming  (17)  to  the  kx,  ky  domain, 


1  f°° 

s(kx,ky)  =  ^  J  h(kx  -  k'T.ky  -  k'y)^(kx,  ky\ k'x,  k'y)r0(k'x.  k'y)  dk'T  dky  , 


:i8) 


where  L  is  defined  bv 
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Definition  of  the  Scattering  Kernel 

For  plane  wave  sources  it  is  possible  to  define  a  factorization  of  the  scattered  field  into  a  product 
of  the  wavenumber  spectrum  of  the  interface  and  a  function  that  is  called  the  scattering  kernel. 
The  scattering  kernel  is  independent  of  the  interface  roughness  function,  and  contains  the  features 
of  the  scattered  field  related  to  the  material  contrast  and  the  source  frequency  and  illumination 
angle.  For  a  plane  wave  source  of  the  form 

lo  (*;,*;)  =  4*  2r  -  kpxxy  -  kp),  (2o) 

equation  (18)  reduces  to 

§.(kx,ky )  =  h(kx  -  kp,ky  -  kp)L{kx,ky;kpx,kl)rp{kp,kl).  (21) 

The  scattered  field  source  term  in  this  case  is  separated  into  a  part  associated  with  a  particular 
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interface  roughness  function,  h(kx  -  k%,ky  —  k (J),  and  a  part  associated  with  the  material  contrast 
and  the  source  frequency  and  illumination  angle,  ^(kx,  ky;  kP)r.o{k%,  k? ).  The  latter  part  is 
designated  as  the  scattering  kernel. 

Knowledge  of  the  scattering  kernel  allows  one  to  evaluate  the  scattering  potential  of  a  model 
independent  of  any  particular  interface.  For  example,  the  transmission  scattering  kernels  for  the 
two-dimensional  model  in  Figure  2  are  given  in  Figure  3.  The  source  is  a  normally  incident,  planar 
P  wrave.  The  scattering  kernels  here  have  been  resolved  into  down-going  P  and  S  wave  scattering 
coefficients  using  a  technique  described  in  the  next  section.  Superimposed  on  these  plots  is  the 
Fourier  transform  of  the  interface.  The  scattered  field  is  simply  the  product  of  the  two  curves. 
From  these  plots  it  is  clear  that  an  interface  with  a  smaller  correlation  length,  and  hence  a  broader 
spectrum  in  the  transform  domain,  would  result  in  large  amplitude  cusps  for  large  scattering  angles 
in  P  and  S.  For  the  transmitted  P  wave,  the  scattered  wave  amplitude  increases  for  scattering  angles 
larger  than  the  P  reflection  critical  angle.  An  amplitude  boost  at  large  angles  may  also  be  seen 
in  the  P  and  S  reflection  coefficients  as  it  will  be  shown  later  with  examples.  This  effect  has  also 
been  demonstrated  by  Levander  and  Hill  (1985)  using  finite  difference  methods  and  by  by  Paul  and 
Campillo  (1988)  using  boundary  integral  equation  methods. 

Describing  Reflection  and  Transmission 

The  source  term  given  by  (18)  generates  P  and  S  waves  above  and  below  the  interface.  To  determine 
the  displacement  coefficients  of  these  waves  requires  a  relation  between  the  displacement-stress 
vector  and  the  up-  and  down-going  P,  SV,  and  SH  wave  components,  the  form  of  which  is  given  by 

r  =  F_b  .  (22) 

where  b  =  [P  S  T  P  S  T)1 ,  the  grave  and  acute  symbols  '•  and  •  denote  down-  and  up-going  waves, 
respectively,  and  P,  S,  and  T  are  the  displacement  coefficients  of  P,  SV,  and  SII  waves,  respectively. 
Using  the  reflection  coefficient  sign  conventions  of  Aki  and  Richards  (1980),  F  is  given  by 
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with  K  =  ^kl  +  k 2,  a  =  s/(X  +  2 n)/p,  0  = 

:  vVp,  7  = 

\JiJl  joi1 

-  A'2,  and  i/  = 

T7.  To 

recover  scattered  field  displacements  from  the  source  term  note  that 


s  =  r (2)  -  I(1)  =  F  (2)6(2)  -  F  (1)6(1)  (24) 

and  that  s  generates  no  down-going  waves  in  the  upper  medium  and  no  up-going  waves  in  the 
lower  medium.  Hence, 
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The  wave  displacement  coefficients  generated  by  s  are  given  by 

b=F=~'s  (26) 
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The  inverse  of  £_  exists  for  all  values  of  ( kx,ky )  except  (kx,ky)  =  (0,0).  At  this  point  in  the  k 
plane  the  scattered  waves  are  propagating  normal  to  the  mean  planar  surface,  and  SV  waves  are 
indistinguishable  from  SH  waves.  For  example,  consider  an  S  plane  wave  traveling  in  the  2  direction 
with  particle  motion  in  the  x  direction.  If  the  wave  direction  is  slightly  perturbed  in  the  x  direction 
it  becomes  an  SV  wave.  A  perturbation  in  the  y  direction  makes  it  an  SH  wave.  These  distinctions 
are  true  even  when  the  perturbations  are  infinitesimal.  Hence,  to  remove  the  singularity  of  £,  an 
arbitrary  naming  convention  must  be  adopted  for  vertically  propagating  S  waves.  In  this  paper, 
the  x  component  of  such  waves  is  labeled  SV,  and  the  y  component  is  labeled  SH.  The  F_  matrix 
for  vertical  propagation  is 
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COMPARISON  WITH  FINITE  DIFFERENCE 


The  simple  form  of  the  scattered  field  source  term  (18)  was  made  possible  by  the  truncation  of 
(15)  to  yield  the  second-order  perturbation  approximation  given  in  (16).  The  error  resulting  from 
the  exclusion  of  higher  order  terms  is  difficult  to  evaluate  analytically.  It  is  possible,  though, 
to  determine  bounds  on  the  domain  of  validity  of  this  approximation.  Kennett  (1972)  derived 
two  conditions  on  the  model  which  must  hold  in  order  for  (18)  to  be  valid.  The  first  condition 
constrains  the  scattered  field  to  be  much  weaker  than  the  background  field,  a  requirement  for  the 
single  scattering  approximation  to  apply.  This  condition  is  expressed  by  the  relation 

s{kx,ky)  <  max|r0(*r,fry)|,  (28) 

fcjAy 

where  s  is  the  scattered  field  source  term  defined  by  (18).  Kennett  (1972)  reduced  this  to  a  simpler, 
but  stricter,  form  by  replacing  the  convolution  integral  in  (18)  by  an  upper-bound  approximation, 
yielding 


KqU!  L 
7T  3] 


r]\ 2  max \h(x)\  <  1, 


max  |/).x(x)|  <  1, 


(29) 
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where 

&\pi  -  OC2P2  PiPi  ~  P2P2 
O  1  p\  +  &2P2  01  Pi  +  02P2 

L  is  the  periodicity  length  of  the  rough  interface,  and  the  wavenumber  spectrum  of  the  incident 
field  is  bounded  by  |/t|  <  k0-  The  second  condition  is  that  the  background  field  must  not  contain 
wavenumbers  so  close  to  grazing  incidence  that  shadow  zones  form.  Shadow  zones  will  be  avoided 
if  the  radius  of  curvature  of  the  interface  is  much  longer  than  a  wavelength.  Such  waves  will  then 
propagate  as  guided  waves  along  the  interface.  This  condition  is  expressed  by 

1  —  max  |h  rjj  <  1,  (30) 

v/«S  +  *?  ' 

where  kz  is  the  vertical  wavenumber  component  associated  with  the  maximum  horizontal  compo¬ 
nent  Kq. 

These  conditions  are  not  satisfactory  for  practical  use,  however.  Approximations  used  in  the 
derivation  of  (29)  result  in  a  much  stricter  bound  than  is  necessary,  and  the  practical  limit  imposed 
by  (30)  needs  to  be  better  defined.  In  order  to  empirically  construct  more  realistic  constraints 
on  interface  roughness,  reflection  and  transmission  coefficients  obtained  using  the  perturbation 
method  described  above  will  be  compared  with  those  derived  from  two-dimensional  finite  difference 
solutions.  By  comparing  results  for  a  range  of  interface  height  and  slope  statistics,  the  domain  of 
validity  of  the  perturbation  method  can  be  explored. 

The  finite  difference  algorithm  used  here  is  a  two-dimensional,  second-order  formulation.  To 
summarize,  the  wave  equations  (1)  and  (2)  are  solved  in  the  time-space  domain  by  replacing  the 
time  and  space  derivatives  by  their  second-order  centered-difference  approximations.  The  accuracy 
is  improved  by  using  a  staggered  mesh  formulation  in  which  horizontal  and  vertical  displacements 
are  represented  on  separate  grids,  each  shifted  by  half  of  the  grid  point  spacing  in  both  coordinate 
directions  with  respect  to  the  other.  This  has  the  added  benefit  of  increasing  the  grid-point  density 
by  a  factor  of  two.  The  formulation  differs  from  that  of  Virieux  (1986)  in  that  we  use  the  second- 
order  displacement/stress  wave  equations  instead  of  the  first-order  velocity/stress  equations.  This 
modification  improves  efficiency,  since  the  final  solution  was  desired  in  terms  of  displacement.  The 
stability  conditions  and  grid  dispersion  relations  are  identical  with  those  of  Virieux’s  formulation. 
A  low-order,  staggered- mesh  scheme  was  chosen  because  it  is  the  most  efficient  method  when  dense 
grid  point  spacings  are  necessary,  as  is  the  case  in  our  models  where  small  interface  irregularities 
must  be  represented.  Higher  order  schemes,  such  as  those  which  use  fourth-order  (Bayliss,  et  0/., 
1980)  or  Fourier  spatial  derivative  operators  (Kosloff  et  al .,  1984;  Fornberg,  1988),  are  generally 
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thought  to  be  more  efficient  than  second-order  schemes,  but  this  is  not  necessarily  true.  The  dense 
grid  point  spacing  used  in  modeling  small  interface  perturbations  puts  the  second-order  derivative 
operator  well  within  its  domain  of  acceptable  accuracy,  and  the  short  operator  length  makes  it  the 
most  efficient  scheme. 

The  finite  difference  method  is  known  to  be  accurate  when  the  wave  fields  and  the  model  are 
both  well  discretized.  Hence,  when  models  with  large  interface  irregularities  are  compared,  the 
finite  difference  method  will  be  used  as  the  standard.  On  the  other  hand,  the  perturbation  method 
is  accurate  for  models  with  small  interface  height  and  slope  irregularities  and  will  be  used  as  the 
standard  for  such  models.  The  comparison  with  finite  difference  solutions  will  proceed  in  two  parts. 
First  it  is  necessary  to  show  that  the  range  of  interface  irregularities  over  which  the  finite  difference 
method  is  accurate  extends  into  the  range  over  which  the  perturbation  method  is  accurate.  This 
will  be  done  by  showing  that  the  finite  difference  method  results  match  the  perturbation  method 
results  for  a  model  with  very  small  height  and  slope  irregularities.  If  the  finite  difference  method 
works  well  in  this  case,  results  for  larger  interface  irregularities  will  also  be  valid  because  they  will 
be  more  accurately  represented  on  the  finite  difference  grid.  Finite  difference  solutions  will  then  be 
used  a  basis  for  comparison  with  perturbation  solutions  for  a  series  of  models  with  larger  height 
and  slope  irregularities  in  order  to  probe  the  limits  of  validity  of  the  perturbation  approximation. 

Accuracy  of  the  Finite  Difference  Method 

All  results  for  scattering  from  a  rough  interface  in  this  paper  are  expressed  as  reflection  and  trans¬ 
mission  coefficients.  The  procedure  for  deriving  these  coefficients  from  the  finite  difference  solution 
is  similar  to  the  method  used  for  deriving  them  for  the  perturbation  solution.  In  short,  equation 
(22)  is  used  to  convert  r  into  b,  and  then  P  and  S  wave  amplitudes  in  b  are  scaled  by  the  source 
wave  amplitude.  The  displacements  and  stresses  in  r  are  computed  by  the  finite  difference  method 
in  the  time-space  domain  along  a  horizontal  linear  array  of  uniformly  spaced  receivers  that  do  not 
intersect  the  interface  at  any  point.  The  receiver  array  should  be  located  between  the  source  and 
the  interface  in  order  to  allow  incident  waves  to  be  separated  from  reflected  waves  based  on  whether 
they  are  traveling  downward  or  upward,  r  is  then  Fourier  transformed  into  the  w-fc  domain  to  yield 
a  form  suitable  for  use  in  (22).  The  distance  of  the  receiver  array  from  the  interface  is  accounted 
for  by  the  z  term  in  the  definition  of  F  in  equation  (23).  If  the  finite  difference  model  is  such 
that  reflections  from  the  top  or  bottom  edges  of  the  grid  will  arrive  within  the  seismogram  time 
window,  absorbing  boundaries  must  be  implemented  which  attenuate  those  reflections  to  a  level 
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much  smaller  than  the  reflection  and  transmission  coefficients  to  he  measured.  Periodic  boundary 
conditions  were  applied  on  the  horizontal  boundaries  of  the  grid  in  order  to  duplicate  the  periodicity 
of  the  spatial  Fourier  transform,  and  the  length  of  the  array,  was  set  to  the  horizontal  dimension 
of  the  finite  difference  grid  to  take  advantage  of  this  periodicity.  For  a  point  source,  L  controls  the 
range  of  incidence  and  scattering  angles  that  are  within  the  receiver  array,  and  the  horizontal  grid 
size  should  be  large  enough  to  capture  the  incidence  and  scattering  angles  of  interest.  Since  the 
range  of  incident  and  scattering  angles  detected  is  also  dependent  on  the  distance  of  the  source  and 
the  receiver  array  from  the  interface,  it  is  best  to  put  the  source  and  the  receiver  array  as  close  to 
the  interface  as  possible  in  order  to  minimize  the  grid  size  and  the  seismogram  length.  Receiver 
spacing  controls  the  maximum  representable  wavenumber  in  r ,  and  was  set  equal  to  the  grid  point 
spacing  to  avoid  aliasing. 

The  reflection  and  transmission  coefficients  obtained  from  finite  difference  modeling  are  com¬ 
pared  with  analytical  coefficients  (Cerveny  ct  al. ,  1977)  for  the  case  of  a  planar  interface  model. 
The  model  is  shown  in  Figure  4.  The  source  has  an  18  Hz  Ricker  time  function  of  the  form 

R(t)  =  (l-  (31) 

where  u>0  is  the  primary  angular  frequency  of  the  wavelet.  The  source  is  implemented  as  a  body 
force  representing  a  point  explosion  source  smoothed  by  a  Gaussian  with  standard  deviation  equal 
to  the  grid  point  spacing.  This  smoothing  is  necessary  to  make  the  infinite  bandwidth  explosion 
source  dipoles  representable  on  the  finite  bandwidth  finite  difference  grid.  This  source  isotropically 
radiates  pure  I’  waves.  It  is  located  20  grid  points  above  the  interface,  and  the  two  arrays  of 
receivers  are  located  10  grid  points  above  and  10  grid  points  below  the  interface.  Reflection  and 
transmission  coefficients  derived  from  the  finite  difference  seismograms  for  this  model  are  shown 
in  Figure  5  along  with  analytical  coefficients.  There  is  generally  good  agreement  for  both  pre-  and 
post-critical  waves.  The  small  disagreement  present  at  the  larger  scattering  angles  results  from  the 
finite  aperture  of  the  receiver  array. 

The  finite  difference  and  perturbation  methods  will  now  be  compared  fora  model  with  very  small 
height  and  slope  perturbations  shown  in  Figure  2.  The  source  is  a  downward  propagating,  planar  P 
wave  with  ar.  18  Hz  Ricker  time  function.  The  interface  roughness  has  a  Gaussian  autocorrelation 
function  with  a  correlation  length  of  l,  =  100  m  and  an  RMS  height  deviation  of  a  =  1(5.7  m.  The 
interface  is  periodic  with  a  period  of  2.70  km,  the  width  of  the  model.  The  finite  difference  sampling 
parameters,  Ax  =  6.(57  m,  A z  =  2.22  m,and  A(  =  O.OOlfifi  s,  result  in  a  maximum  phase  dispersion 
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error  of  —1.1  percent  at  18  Hz.  The  two  receiver  arrays  are  located  as  close  to  the  interface  as 
is  possible  without  intersection.  The  Fourier  transform  of  the  interface  function  which  has  been 
discretized  for  the  finite  difference  grid  is  shown  in  Figure  6a  as  a  function  of  horizontal  slowness 
evaluated  at  18  Hz.  The  zero  spectral  power  at  the  origin  results  from  the  zero  mean  interface 
deviation.  Histograms  of  interface  height  and  slope  are  shown  in  Figures  6b  and  6c.  The  largest 
deviation  from  the  mean  planar  surface  is  12  m,  or  15  percent  of  the  S  wavelength  in  the  upper 
medium.  The  largest  slope  is  0.16,  while  the  mean  and  standard  deviation  of  the  slope  are  zero  and 
0.059,  respectively.  Since  the  maximum  interface  height  and  slope  are  fairly  well  approximated  by 
twice  their  standard  deviations,  histograms  will  not  be  provided  for  the  remaining  rough  interfaces 
used  in  this  study.  A  comparison  of  the  reflection  and  transmission  coefficients  derived  from  the 
finite  difference  and  perturbation  methods  is  shown  in  Figure  7.  Agreement  is  very  good,  with 
both  amplitude  and  shape  predicted  well  by  the  perturbation  method.  “Ringing”  appears  to  some 
degree  on  all  of  the  finite  difference  coefficient  plots,  and  is  most  apparent  on  the  S  wave  coefficient 
plots.  Another  check  on  the  accuracy  of  the  reflection  and  transmission  coefficient  results  is  to  plot 
the  coefficients  for  the  non-physical  waves  in  the  solution:  the  down-going  S  in  the  upper  medium 
and  the  up-going  S  and  P  in  the  lower  medium.  Plots  of  the  non-physical  coefficients  for  the  finite 
difference  data  are  shown  in  Figure  8.  In  general,  the  amplitude  of  the  non-physical  waves  increases 
with  increasing  scattering  angle,  since  the  wavenumbers  corresponding  to  these  larger  angles  are 
less  well  represented  by  the  receiver  array.  The  down-going  P  wave  in  the  upper  medium  has  unit 
amplitude  at  zero  angle.  It  is  clear  from  this  figure  that  the  maximum  error  is  less  than  0.1  percent, 
and  that  this  error  is  associated  with  the  S  wave  in  the  upper  medium.  This  error  is  enough  to 
explain  the  difference  between  the  finite  difference  and  perturbation  results  for  the  reflected  S  wave 
coefficient  at  large  scattering  angles. 

Accuracy  of  the  Perturbation  Method  in  the  u>-k  Domain 

The  domain  of  validity  of  the  perturbation  method  will  be  explored  by  comparing  reflection  and 
transmission  coefficients  generated  using  the  perturbation  method  with  those  derived  from  finite 
difference  modeling.  Six  rough  interface  models  are  used  in  this  comparison,  each  having  a  Gaus¬ 
sian  autocorrelation  function  and  a  uniform  random  phase.  The  interface  functions  are  shown  in 
Figure  9(a-f)  with  two  times  vertical  exaggeration.  RMS  interface  slope  ranges  from  0.037  to  0.99 
and  RMS  interface  height  is  0.01  km  for  models  A-E  and  0.015  Km  for  model  F.  Since  the  accuracy 
of  the  perturbation  method  is  sensitive  to  the  smallest  elastic  wavelength  in  the  model,  interface 
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height  is  measured  here  in  terms  of  S  wavelengths  in  the  upper  medium  (Si),  and  thus  the  relative 
height  of  the  irregularities  changes  as  the  frequency  changes.  The  bandwidth  of  the  18  Hz  Ricker 
wavelet  time  function  used  in  the  finite  difference  calculation  allowed  reflection  and  transmission 
coefficients  to  be  computed  at  three  frequencies  (9.93,  16.5,  and  26.5  Hz)  in  each  of  the  six  models, 
for  a  total  of  18  comparisons  for  each  coefficient.  In  these  comparisons,  RMS  interface  height  ranges 
from  0.069  to  0.28  Si  wavelengths.  RMS  height,  correlation  length,  and  RMS  slope  for  the  models 
used  are  listed  in  Table  2.  The  material  parameters  are  Qi  =  2.50  km/s,  0i  -  1.44  km/s,  py  =  1.00 
g/cm3,  a2  =  3.00  km/s,  /32  =  1.73  km/s,  and  p2  =  1.00  g/cm3,  the  same  as  in  Figure  2.  The  source 
is  a  normally  incident  planar  P  wave  in  layer  one.  Plots  of  t  ese  comparisons  are  shown  for  the 
selected  models  A,  D,  and  E  in  Figures  lO(a-c).  The  reflection  and  transmission  coefficients  are 
the  displacement  coefficients  Pi,  $i,  P2,  and  52  are  normalized  so  that  P(  =  1, 

In  order  to  facilitate  the  comparison  of  the  perturbation  and  finite  difference  coefficients,  the 
curves  are  compared  by  finding  the  L2  norm  difference  defined  by 


II  Ppt  ~  Pjdh  =  100  x 


Jo f[Ppt(0,)  -  Pfd(9s)}2dds 

f}  PJd{9,yd6, 


(32) 


where  6S  is  the  scattering  angle  and  the  P  coefficient  is  chosen  for  illustration.  The  differences 
between  method  results  are  plotted  for  each  of  the  four  coefficients  in  Figures  ll(a-d)  for  a  constant 
RMS  interface  height  of  0.01  km,  and  in  Figures  ll(e-h)  for  a  constant  RMS  interface  slope  of  0.1. 
I  he  constant  height  plots  show  that  the  accuracy  of  the  perturbation  solution  is  not  significantly 
degraded  as  RMS  interface  slope  varies  within  the  range  tested.  This  L2  norm  result  is  verified 
in  the  coefficient  comparison  plots  in  Figures  10(a-d).  The  constant  slope  plots  show  that  error 
increases  rapidly  with  increasing  RMS  interface  height,  and  in  this  particular  example  is  acceptable 
for  values  of  RMS  height  less  than  about  0.1  Si  wavelengths. 


Accuracy  of  the  Perturbation  Method  in  the  t-x  Domain 

Error  in  time  domain  seismograms  generated  using  the  perturbation  approximation  cannot  be 
directly  estimated  from  the  error  results  presented  for  the  u i-k  domain.  This  is  because  each 
seismogram  contains  energy  from  a  broad  spectrum  of  frequencies  and  wavenumbers.  Here  the  time 
domain  error  will  be  estimated  by  comparing  reflection  seismograms  generated  by  the  perturbation 
method  with  those  generated  from  the  finite  difference  method.  To  generate  reflection  seismograms 
using  the  perturbation  method,  the  scattering  coefficients  evaluated  using  the  two-dimensional  form 
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of  (26)  are  converted  to  displacement  by  using  (22)  with  ky  —  0.  These  displacements  are  then 
Fourier  transformed  in  u  and  k  to  yield 


ux 

uz 


1  r°° 

4r2  J-oo 
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e,kxdk 

-OO 


F\4  *15 
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(33) 


Fourier  transforms  are  evaluated  using  the  discrete  wavenumber  method  (Bouchon,  1977).  Seismo¬ 
grams  generated  with  (33)  contain  only  the  scattered  field.  The  background  field  can  be  added  if 
desired. 

Seismograms  were  generated  for  each  of  the  interface  functions  parameterized  in  Table  3  using 
both  the  perturbation  and  finite  difference  methods,  and  a  subset  of  these  are  shown  overlain  in 
Figure  12.  The  source  and  receiver  configurations  and  the  material  parameters  are  the  same  as  those 
used  in  the  last  section.  The  source  time  function  is  an  18  Hz  Ricker  wavelet.  These  seismograms 
show  that  the  shapes  of  the  waveforms  generated  by  the  perturbation  and  finite  difference  methods 
agree  fairly  well  for  the  range  of  models  considered.  Errors  in  amplitude  and  traveltime  seem  to  be 
positively  correlated,  both  being  larger  at  receivers  that  are  over  extrema  in  the  interface  height 
function  than  at  other  receiver  locations.  This  feature  will  be  explored  next  using  an  L2  norm  error 
measure.  Since  the  total  waveform  is  composed  of  the  background  field  added  to  the  scattered  field, 
total  field  perturbation  method  seismograms  minus  the  source  waves  are  presented  for  comparison 
in  Figure  13  for  models  A-I. 

The  L2  norm  was  calculated  as  a  measure  of  error  between  the  perturbation  and  finite  difference 
seismograms,  using 


E(x,At)  =  ||uz(x)pi  -  uJ(x)/c'|j2(At) 


=  100  x 


(34) 


where  L  is  the  length  of  the  receiver  array  and  At  is  the  time  shift  between  the  two  seismograms. 
The  numerator  was  efficiently  evaluated  using  the  fast  Fourier  transform  since  the  seismograms 
are  periodic  in  space.  The  denominator  is  the  average  power  of  the  entire  array  of  finite  difference 
seismograms  in  order  to  minimize  the  effect  of  low  amplitude  seismograms  on  E(At).  If  E(At) 
is  normalized  instead  to  the  current  seismogram,  globally  small  error  can  overwhelm  such  low 
amplitude  seismograms  to  produce  large  values  of  E(At). 
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E(At )  is  used  to  provide  two  separate  measures  of  error  for  each  seismogram:  travel  time 
error,  which  is  defined  as  the  time  delay  At  at  which  E  is  minimum;  and  amplitude  error,  which 
is  defined  as  E  =  E(At).  A  plot  of  At  versus  seismogram  offset  x  for  Model  A  is  shown  in 
Figure  14.  Overlain  on  this  plot  is  a  plot  of  the  one-way,  vertical  P  wave  travel  time  associated 
wdth  the  interface  function  In  the  absence  of  travel  time  error,  At  would  be  zero  for  all  the 

seismograms.  The  agreement  between  these  two  curves  indicates  that  travel  time  from  the  interface 
to  the  receiver  in  the  perturbation  method  is  referenced  to  the  mean  planar  interface,  and  not  to 
the  perturbed  interface. 

The  reason  for  the  traveltime  error  can  easily  be  seen  for  the  case  where  the  interface  perturba¬ 
tion  function  h{x)  is  a  constant.  In  this  case,  the  perturbation  expansion  given  by  (15)  transforms 
to  an  exact  form  in  the  u-k  domain: 


r(u,kx,ky;h)  =  L{L  +  hA  +  -yA2  +  ...)r_{u,kx,ky\Q)  (35) 

=  eh=r(u>,kx,ky;  0)  (36) 

=  Fek/i£-lr(v,kxtkv;0).  (37) 

Equation  (36)  is  simply  the  definition  of  a  propagator  matrix  which  transforms  the  solution  for  r 
at  depth  zero  to  the  solution  at  depth  h.  Equation  (37)  uses  a  diagonal  factorization  of  A  of  the 
form 


A  =  £A£-1  (38) 

where  £_  is  defined  by  (23)  and  is  always  invertible,  and  A  is  a  diagonal  matrix  whose  diagonal 
elements  in  the  three-dimensional  case  are  An  -  -A44  =  ij  and  A22  =  A33  =  -A55  =  -A66  =  iv- 
The  exact  form  of  the  source  term  (18)  valid  for  constant  h  is  then 

s{u,ki,ky)  =  (eh-  *  -  ek-7))r{u,kx,ky;0)  (39) 

=  (£(1)efc4(,,F(1)-1  -  L{2)eh=WF}2)~l)r(uykx,ky;0)  (40) 

»  (£(1)(/i4(1))£(1)_1 -£(2)(/iA(2))£(2)"l)r(a;,itI,Jby;0)  (41) 

+  0  (h2). 

I  lie  phase  factors  of  the  form  eh=  in  (40)  are  responsible  for  adjusting  the  travel  time  to  account 
for  Uie  height  of  the  interface.  The  accuracy  of  the  first  order  approximation  eh=  —  /  ~  hA  in  (-11) 
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controls  the  accuracy  of  the  travel  time  adjustment  and  influences  the  accuracy  of  the  scattered 
field  amplitude.  The  accuracy  of  the  first  and  second  order  approximations  to  e'h=  —  /  is  illustrated 
in  Figure  15.  This  figure  shows  that  at  least  a  second  order  approximation  is  required  for  accurate 
representation  of  phase,  and  therefore  of  traveltime,  although  both  approximations  converge  to  the 
exact  solution  as  h  — < ►  0.  However,  the  modulus  of  the  first  order  approximation  is  more  accurate 
than  that  of  the  second  order  approximation,  but  this  is  a  consequence  of  the  non-monotonic 
improvement  of  accuracy  as  approximation  order  increases. 

RMS  amplitude  error  E  in  the  comparison  of  perturbation-  and  finite  difference-derived  seis¬ 
mograms  is  shown  in  Figure  16.  RMS  errors  of  the  seismograms  of  each  model  were  averaged,  and 
these  mean  errors  are  plotted  against  RMS  height  for  RMS  slopes  of  0.037  and  0.10,  and  against 
slope  ”  an  RMS  height  of  0.125  S i  wavelengths.  Contrary  to  the  results  obtained  in  the  u i-k 
domain,  RMS  amplitude  error  clearly  increases  with  increasing  RMS  interface  slope,  as  well  as 
with  increasing  RMS  interface  height.  Comparison  of  the  rates  of  error  growth  for  the  fixed  slope 
and  fixed  height  error  plots  shows  that  error  increases  with  increasing  height  nearly  twice  as  fast 
as  it  does  with  increasing  slope. 

When  the  time  shift  error  is  removed,  amplitude  error  is  sensitive  to  two  factors:  the  amplitude 
scale  factor  and  the  waveform  character.  Waveform  character  refers  to  the  general  shape  of  the 
waveform  including  the  presence  and  length  of  the  coda.  It  is  fairly  well  approximated  in  all  of 
the  nine  models  considered  here.  With  the  exception  of  small  amplitude  arrivals,  the  shapes  of  the 
first  arrivals  generally  agree  well  with  the  finite  difference  results,  even  for  models  with  RMS  slopes 
as  high  as  one,  and  the  complexity  and  decay  rate  of  the  coda  are  also  well  approximated.  The 
amplitude  scale  factor,  however,  is  generally  overestimated  by  the  perturbation  method.  This  is 
true  to  the  nature  of  the  Born  approximation,  which  violates  conservation  of  energy  by  introducing 
the  rough  interface  effects  via  a  source  term  whose  energy  is  added  to  the  energy  of  the  background 
field.  With  these  points  in  mind,  the  range  of  validity  spans  the  entire  range  of  the  models  tested 
if  the  error  criterion  is  based  on  waveform  character.  On  the  other  hand,  if  the  error  criterion  is 
based  on  amplitude  scale  or  the  absolute  RMS  error  the  range  is  limited  to  RMS  slopes  of  less 
than  about  0.3  and  RMS  heights  of  less  than  about  0.25  of  the  smallest  wavelength  present  in  the 
scattered  field.  Notice  that  the  latter  criterion  results  is  a  less  restrictive  range  than  the  RMS  error 
analysis  on  scattering  coefficients  done  in  the  u-k  domain. 
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FEATURES  OF  THE  3-D  SCATTERED  FIELD 


Three-dimensional  scattering  examples  will  now  be  presented  for  the  model  shown  in  Figure  17. 
The  material  and  interface  parameters  are  identical  to  those  of  model  D  (see  Table  2  and  Figure  9) 
with  the  Gaussian  auto-correlation  function  having  identical  correlation  lengths  in  the  x  and  y 
directions.  In  general,  a  Gaussian  rough  interface  whose  major  and  minor  axes  coincide  with  the  x 
and  y  axes  has  a  Gaussian  auto-correlation  function  of  the  form 


{kx^ky^j 


~  e  < 


(42) 


where  Lx  and  Ly  are  the  major-  and  minor-axis  correlation  lengths.  Equation  (42)  can  be  gener¬ 
alized  to  allow  roughness  trends  at  an  angle  9  to  the  x  axis  by  applying  a  rotation  transformation 
to  get 


h(kx,ky)h*(kx,ky)  ~  e- »!„(«)+<:„  cos(fl))2]i  (43) 

In  the  first  example,  the  source  is  an  18  Hz  planar  SV  wave  propagating  downward  along  the  2  axis 
with  particle  motion  in  the  x  direction.  The  SV  wave  scatters  into  reflected  and  transmitted  P, 
SV,  and  SH  waves.  The  three-dimensional  P,  SV,  and  S II  transmission  coefficients  for  this  model, 
found  using  the  perturbation  method,  are  given  in  Figure  18.  The  x  and  y  scattering  angles  in 
these  figures,  referred  to  as  4>x  and  <py,  are  measured  from  the  downward  2  axis  in  the  x-z  and 
y-z  planes,  respectively.  They  are  defined  by  <j)x  =  sin-1  ( kxv/u> )  and  <j>y  =  sin-1  (kyv/u)  where 
v  is  the  body  wave  velocity  of  the  transmitted  wave  concerned.  Scattering  in  the  2-direction,  for 
example,  is  given  by  <px  =  <f>y  =  0.  The  transmission  coefficient  plots  cover  scattering  angles  in  the 
range  -9U°  <  4>x  <  90°  and  0°  <  <f>y  <  90°,  and  are  symmetric  about  the  <j>y  -  0  axis.  For  all  three 
scattered  wavetypes,  these  plots  show  that  scattering  is  maximal  in  the  direction  that  conserves 
source  particle  motion:  for  an  SV  source  with  particle  motion  in  the  x-z  plane,  P  and  SV  are 
maximally  scattered  in  the  x-z  plane,  and  SH  is  maximally  scattered  in  the  y-z  plane.  This  effect  is 
exaggerated  by  the  use  of  the  single  scattering  approximation.  A  null  is  present  in  the  plane  normal 
to  maximal  plane:  the  y-z  plane  for  P  and  SV  waves  and  the  x-z  plane  for  SH  waves.  An  alternative 
display  of  the  three-dimensional  transmission  coefficients  is  a  cross-section  of  the  coefficients  for  all 
scattering  angles  at  a  particular  azimuth,  where  azimuth  is  measured  clockwise  in  the  x-y  plane 
from  the  positive  z-axis,  and  the  scattering  angle  is  defined  by  i/>  =  sin-1  (kzv/u).  Cross-sections 
of  the  reflection  and  transmission  coefficients  for  the  same  model  and  source  parameters  are  given 
for  azimuthal  angles  of  10°,  45°,  and  80°  in  Figure  19. 
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Transmission  scattering  kernels  for  the  model  in  Figure  17  are  given  in  Figure  20.  Comparison 
with  the  SY  and  SH  transmission  coefficients  in  Figure  18  shows  that  the  spatially  band-limited 
nature  of  this  particular  interface  damps  out  the  large  amplitude  features  present  at  large  scat¬ 
tering  angles.  The  P  transmission  coefficient  exhibits  a  hint  of  the  cusps  present  in  the  kernel. 
Consideration  of  another  interface  roughness  function  requires  only  a  visual  superposition  of  the 
new  interface  spectrum  with  the  kernel.  For  the  same  SV  source  as  above  and  a  two-dimensional 
rough  interface  with  variation  in  the  x  direction,  the  problem  is  truly  two-dimensional,  and  waves 
are  scattered  into  P  and  SV.  If  the  interface  is  rotated  90  degrees  so  that  variation  is  in  the  y 
direction,  the  problem  is  fully  three-dimensional,  and  waves  are  scattered  into  P  and  SH.  More 
general  interfaces  whose  auto-correlation  functions  are  described  by  (43),  or  perhaps  a  von  Karman 
or  exponential  function,  are  handled  with  the  same  approach. 

The  second  example  is  the  same  as  the  first,  but  with  the  SV  source  replaced  by  a  P  source.  The 
three-dimensional  P  and  SV  transmission  coefficients  for  this  model  are  given  in  Figure  21.  Since 
the  source  particle  motion  and  the  Gaussian  auto-correlation  function  are  azimuthally  invariant, 
the  scattering  coefficients  are  also  azimuthally  invariant.  Deviations  of  the  contours  from  circular 
arcs  are  artifacts  of  the  contouring  program.  Reflection  and  transmission  coefficient  cross-sections 
for  this  model  are  given  in  Figure  22.  Note  that  scattered  SH  waves  are  not  generated  in  this 
example.  This  is  a  consequence  of  the  single-scattering  approximation.  For  a  normal  incidence 
planar  P  wave,  a  minimum  of  two  bounces  are  required  to  generate  SH. 

DISCUSSION  AND  CONCLUSIONS 

A  perturbation  method  has  been  presented  here  for  computing  three-dimensional  body  wave  scat¬ 
tering  from  a  rough  interface.  Its  speed  and  simplicity  are  such  that  each  of  the  examples  in 
this  paper  was  generated  on  a  Macintosh  SE  computer  in  less  than  two  minutes.  Speed  is  an 
important  consideration  when  three-dimensional  modeling  is  necessary.  For  example,  many  of  the 
two-dimensional  finite  difference  computations  done  for  comparison  required  15  Mbytes  of  core  and 
23  hours  of  CPU  time  on  a  Vax  8800.  In  three  dimensions,  finite  difference  solutions  for  these 
models  are  beyond  our  resources  at  present.  For  the  class  of  irregular  interface  models  with  small 
RMS  height  deviations,  the  perturbation  method  is  a  useful  alternative. 

Since  interface  height  deviations  in  the  models  presented  here  are  small,  comparisons  of  the 
perturbation  method  with  the  finite  difference  method  were  preceded  by  careful  testing  of  the  finite 
difference  method  to  show  that  it  is  valid  for  the  small  deviations  used  in  these  comparisons.  It  can 
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be  shown  that  in  the  limit  as  the  finite  difference  grid  sampling  interval  goes  to  zero,  and  as  the 
numerical  precision  of  the  computer  goes  to  infinity,  the  finite  difference  solution  converges  to  the 
exact  solution  as  O(Ai)  (Brown,  1984).  However,  in  order  for  the  grid  sampling  interval  to  be  small 
enough  that  the  finite  difference  solution  is  sufficiently  close  to  the  exact  solution,  the  number  of 
grid  points  describing  a  fixed  model  must  increase  as  the  size  of  the  interface  height  perturbations 
decreases.  Therefore,  the  core  size  and  speed  of  a  computer  are  constraints  in  the  minimum  interface 
height  perturbation  that  can  be  accurately  modeled  using  finite  difference  method.  Of  course, 
a  planar  boundary  is  an  exception  to  this  limit.  Another  lower  limit  on  interface  perturbation 
size  is  imposed  by  numerical  precision.  As  the  size  of  interface  height  perturbations  decreases, 
the  amplitude  of  scattered  seismic  waves  decreases.  Since  these  scattered  waves  are  added  to 
the  relatively  large  planar  interface  response,  insignificant  numerical  noise  in  the  planar  interface 
response  can  be  significant  when  compared  with  the  scattered  field.  This  is  probably  the  cause  of 
the  “ringing”  present  in  the  finite  difference-derived  scattering  coefficients.  Comparisons  of  finite 
difference-derived  scattering  coefficients  for  a  planar  model  against  an  analytical  solution,  and  for  a 
model  with  small  interface  height  and  slope  perturbations  against  a  perturbation  method  solution, 
show  that  this  “ringing”  tends  to  oscillate  about  the  exact  solution.  Hence,  the  smoothed  finite 
difference  scattering  coefficients  are  useful  for  comparisons. 

Numerous  two-dimensional  model  comparisons  of  finite  difference  and  perturbation  method 
scattering  coefficients  were  made  in  order  to  determine  the  range  of  validity  of  the  perturbation 
method.  The  L2  norms  of  the  differences  were  computed  for  varying  RMS  slope  and  constant  RMS 
height,  and  for  varying  RMS  height  and  constant  RMS  slope.  There  are  two  major  trends  in  the 
L2  norm  results.  Firstly,  error  increases  strongly  with  increasing  RMS  height,  with  acceptable 
levels  for  RMS  heights  of  less  than  about  0.1  shear  wavelengths.  Secondly,  error  appears  to  be 
roughly  constant  for  increasing  RMS  slope  for  the  tested  range  of  from  0.037  to  0.99.  Smaller 
trends,  such  as  the  apparent  increase  in  accuracy  with  increasing  RMS  slope,  are  misleading.  The 
L2  norm  is  overly  sensitive  to  the  noise  in  the  finite  difference  solution  being  used  as  the  standard 
for  comparison.  This  is  apparent  when  the  scattering  coefficient  plots  are  examined  for  the  trend 
seen  in  the  L2  norm  plots.  The  Li  norm  was  also  calculated  to  see  whether  it  is  loss  sensitive  to 
this  noise,  but  the  improvements  were  minimal.  Ultimately,  all  trends  must  be  confirmed  by  the 
scattering  coefficient  plots. 

Time  domain  seismograms  generated  by  the  perturbation  method  were  compared  with  finite 
difference  seismograms  for  the  same  models.  The  perturbation  method  was  able  to  accurately 
predict  waveform  character  for  the  entire  range  of  models  considered,  including  the  shape  of  the 
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first  arrival  and  the  presence,  duration,  and  complexity  of  the  coda.  Using  RMS  error  analysis 
to  determine  traveltime  error,  it  was  shown  that  when  first  order  perturbation  theory  is  used,  the 
travel  time  from  the  interface  to  the  receiver  is  referenced  to  the  mean  planar  interface,  and  not 
from  the  perturbed  interface.  Hence,  traveltime  error  increases  linearly  with  the  height  of  the 
interface.  It  was  shown  that  a  second  order  theory  is  required  for  travel  time  to  include  the  height 
of  the  perturbation.  With  this  travel  time  error  removed,  the  mean  RMS  error  was  determined  for 
the  perturbation  seismograms  in  each  of  the  models.  If  this  mean  RMS  error  is  the  criterion  for 
determining  the  domain  of  validity  of  the  perturbation  method,  it  is  valid  for  RMS  slopes  of  less 
than  about  0.25  and  RMS  heights  of  less  than  about  20  percent  of  the  smallest  wavelength  in  the 
scattered  field. 

Three-dimensional  scattering  kernels  generated  for  an  SV  plane  wave  normally  incident  on  a 
rough  interface  shows  that  waves  are  maximally  scattered  in  directions  for  which  the  scattered 
wave  particle  motion  coincides  with  that  of  the  incident  wave.  A  P  wave  in  the  same  geometry 
induces  aziinuthally  isotropic  radiation.  In  the  three-dimensional  scattering  example  it  was  shown 
that  there  is  a  null  in  P  and  SV  scattering  in  the  y  direction,  and  a  null  in  SH  scattering  in  the 
x  direction.  These  nulls  do  not  exist  in  the  time-space  domain,  where  a  receiver  in  any  location 
can  detect  waves  traveling  in  all  directions.  These  scattering  kernels  also  show  that  scattered  wave 
amplitudes  tend  to  increase  for  scattering  angles  beyond  the  P  wave  critical  angle  (defined  with 
respect  to  the  planar  interface  model).  This  confirms  the  result  of  Paul  and  Campillo,  1988. 
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Symbol 

Definition 

Equation 

Introduced 

Pj 

Density  in  layer  j 

(1) 

ua 

Angular  frequency 

(1) 

Uj 

j  tli  component  of  displacement 

(1) 

r 

Cauchy  stress  tensor 

(1) 

i  Pi 

Lame  parameters  for  layer  j 

(2) 

r 

Displacement-stress  vector 

(3) 

A 

Wave  equation  coefficient  matrix  for  layer  j 

(3) 

Oj 

Partial  derivative  vv.r.t.  jth  component 

(4) 

0 

Simplification  variable  for  layer  j 

(4) 

Zifa-,  y) 

Unit  normal  to  interface 

(5) 

h{x,y) 

Zero  mean  interface  height  function 

(5) 

Tj 

jth  component  of  traction 

(6) 

Q  x  \  Q  yJ^  Rotation  matrices  for  layer  j 

(7) 

rU) 

Rotated  displacement-stress  vector  in  layer  j 

(7) 

A' 

&th  order  scattered  field  for  layer  j 

(14) 

La 

Displacement-stress  vector  for  mean  planar  interface 

(14) 

s 

Scattered  field  source  term 

(18) 

4 

Kernel  matrix  for  scattered  field  source  term 

(18) 

ky 

x  and  y  components  of  wavenumber 

(18) 

b 

Wave  coefficient  vector 

(22) 

£ 

Layer  matrix 

(22) 

K 

Magnitude  of  horizontal  wavenumber 

(23) 

Qi i  0] 

Compressional  and  shear  wave  speeds  in  layer  j 

(23) 

'll 

Vertical  component  of  wavenumber  for  a  compressional  wave 

(23) 

vi 

Vertical  component  of  wavenumber  for  a  shear  wave 

(23) 

b 

Scattering  coefficient  vector 

(25) 

t 

Scattering  coefficient  transformation  matrix 

(25) 

Ko 

Maximum  bound  on  source  horizontal  wavenumber 

(29) 

9, 

Scattering  angle 

(32) 

A 

Diagonal  eigenvalue  matrix  for  A 

(38) 

Table  1:  Table  of  symbols  used  in  this  paper. 
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Model 

Correlation 

Length 

(km) 

RMS 

Slope 

RMS 

Height 

(km) 

RMS  Height 

in 

Si  Wavelengths 

at  9.93  Hz 

RMS  Height 

in 

Si  Wavelengths 

at  16.5  Hz 

RMS  Height 

in 

Si  Wavelengths 

at  26.5  Hz 

A 

0.010 

0.30 

0.037 

0.069 

0.11 

0.18 

B 

0.010 

0.10 

0.10 

0.069 

0.11 

0.18 

C 

0.010 

0.050 

0.20 

0.069 

0.11 

0.18 

D 

0.010 

0.033 

0.30 

0.069 

0.11 

0.18 

E 

0.010 

0.010 

0.99 

0.069 

0.11 

0.18 

F 

0.015 

0.15 

0.10 

0.10 

0.17 

0.28 

Table  2:  Interface  parameters  for  the  rough  interfaces  used  in  the  u>-k  domain  comparisons.  The 
interface  height  functions  shown  plotted  in  Figure  9. 


Model 

Correlation 

Length 

(km) 

RMS 

Slope 

RMS 

Height 

(km) 

RMS  Height 

in 

Si  Wavelengths 

at  18  Hz 

A 

0.30 

0.037 

0.010 

B 

0.10 

0.10 

0.010 

C 

0.050 

0.20 

0.010 

0.125 

D 

0.033 

0.30 

0.010 

0.125 

E 

0.010 

0.99 

0.010 

0.125 

F 

0.15 

0.10 

0.015 

0.188 

G 

0.45 

0.038 

0.015 

0.188 

H 

0.60 

0.039 

0.020 

0.25 

I 

0.20 

0.10 

0.020 

0.25 

Table  3:  Interface  parameters  for  the  rough  interfaces  used  in  the  t-x  domain  comparisons. 
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FIGURE  CAPTIONS 


Figure  1:  Geometry  of  the  rough  interface  model.  The  interface  is  defined  by  z  —  h(x,y).  and  the 
downward-pointing  unit  normal  to  this  surface  at  each  point  is  denoted  by  n(x,y).  The  2  =  0 
plane  is  defined  as  the  mean  planar  surface  through  the  rough  interface.  <*i,  pi,  »2i  >h,  and 
/>2  denote  the  compressional  and  shear  wave  speeds  and  density  for  the  materials  above  and 
below  the  interface,  respectively. 

Figure  2:  Two-dimensional  rough  interface  model.  The  source  is  a  normally  incident,  18  Hz,  plane 
wave.  The  rough  interface  has  a  Gaussian  autocorrelation  function  with  a  correlation  length  of 
L  =  100  m  and  an  RMS  height  deviation  of  a  =  16.7  m.  The  interface  is  periodic  with  a  period 
of  2.70  km,  the  width  of  the  model  above. 

Figure  3:  P  and  S  wave  transmission  scattering  kernels  for  a  normally  incident  P  wave  for  the 
model  in  Figure  2.  The  Fourier  transform  of  the  interface,  mapped  into  scattering  angle,  is 
shown  superimposed. 

Figure  4:  Planar  interface  model  used  to  compare  finite  difference  derived  reflection  and  transmis¬ 
sion  coefficients  with  analytical  solutions.  The  source  is  a  point  explosion  with  an  18  Hz  Ricker 
wavelet  time  function.  The  sampling  parameters,  Ax  =  Az  =  6.67  m  and  At  =  0.00156  s, 
result  in  a  maximum  phase  dispersion  error  of  -1.1  percent  at  18  Hz.  The  source  is  located 
20  grid  points  above  the  interface,  and  the  receiver  arrays  are  located  10  grid  points  on  either 
side  of  the  interface.  The  horizontal  receiver  interval  is  two  grid  points.  All  waves  within  the 
seismogram  time  window  are  contained  within  the  finite  difference  grid,  eliminating  the  need 
for  absorbing  boundaries. 

Figure  5:  Comparison  between  P  and  S  wave  reflection  and  transmission  coefficients  generated  by 
finite  difference  and  analytic  methods  for  a  planar  interface.  The  model  is  shown  in  Figure  4. 
The  small  disagreement  present  at  the  larger  scattering  angles  is  “Gibb’s  ringing”  that  results 
from  the  finite  aperture  of  the  receiver  array. 

Figure  6:  Properties  of  of  the  interface  shown  in  Figure  2:  (a)  Fourier  transform,  where  the 
horizontal  slowness  is  evaluated  for  18  Hz,  and  histograms  of  (b)  interface  height  and  (c)  interface 
slope. 

Figure  7:  Comparison  between  P  and  S  wave  reflection  and  transmission  coefficients  generated  from 
the  finite  difference  and  perturbation  methods  for  the  model  shown  in  Figure  2. 

Figure  8:  Amplitudes  of  waves  incident  on  the  interface  corresponding  to  the  scattered  wave  am¬ 
plitudes  shown  in  Figure  7.  These  are  provided  as  a  measure  of  the  error  versus  angle.  The 
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down-going  P  wave  in  the  upper  medium  has  unit  amplitude  at  zero  angle. 

Figure  9:  Interface  functions  used  in  comparison  of  reflection  and  transmission  coefficients  derived 
from  the  perturbation  method  with  those  derived  from  the  finite  difference  method.  The  inter¬ 
faces  have  a  Gaussian  autocorrelation  function  with  RMS  height,  correlation  length,  and  RMS 
slope  for  each  interface  listed  in  Table  2.  The  functions  are  displayed  with  two  times  vertical 
exaggeration. 

Figure  10a:  Comparison  of  reflection  and  transmission  coefficients  derived  from  the  perturbation 
and  finite  difference  methods  for  model  A.  The  parameters  of  the  rough  interface  are  given  in 
Table  2  and  the  interface  function  is  illustrated  in  Figure  9. 

Figure  10b:  Comparison  of  reflection  and  transmission  coefficients  derived  from  the  perturbation 
and  finite  difference  methods  for  model  D.  The  parameters  of  the  rough  interface  are  given  in 
Table  2  and  the  interface  function  is  illustrated  in  Figure  9. 

Figure  10c:  Comparison  of  reflection  and  transmission  coefficients  derived  from  the  perturbation 
and  finite  difference  methods  for  model  E.  The  parameters  of  the  rough  interface  are  given  in 
Table  2  and  the  interface  function  is  illustrated  in  Figure  9. 

Figure  11:  L2  norm  comparison  of  finite  difference  and  perturbation  method  derived  reflection  and 
transmission  coefficients,  (a-d)  RMS  interface  height  is  a  constant  0.01  km  and  RMS  slope 
varies  from  0.037  to  0.99.  RMS  interface  height  for  three  frequencies  can  be  expressed  as  0.069, 
0.11,  and  0.18  Si  wavelengths,  (e-h)  RMS  interface  slope  is  a  constant  0.1  and  with  RMS 
interface  height  varies  from  0.069  to  0.28  Si  wavelengths. 

Figure  12:  Representative  scattered  field  seismograms  generated  for  each  of  the  interface  functions 
parameterized  in  Table  3.  The  source  is  a  normally  incident  plane  wave  with  an  18  Hz  Ricker 
time  function.  The  bold  curves  are  finite  difference  solutions,  and  the  lighter  curves  are  pertur¬ 
bation  solutions.  The  four  vertical  component  receivers  are  0.0666  km  above  the  mean  planar 
interface  at  horizontal  offsets  of  0.0133,  1.35,  2.68,  and  4.02  km.  All  seismograms  are  plotted 
at  the  same  scale. 

Figure  13:  The  seismograms  of  Figure  12  that  were  generated  by  the  perturbation  method  are 
shown  here  with  the  background  reflected  field  included  in  order  to  show  the  total  waveform 
(minus  the  source  wave).  Plotted  at  the  same  scale  as  Figure  12. 

Figure  14:  Travel  time  error  At  plotted  against  seismogram  offset  x  for  Model  A.  Shown  overlain  is 
a  plot  of  the  one-way,  vertical  P  wave  travel  time  associated  with  the  interface  function  <*?)■ 

Figure  15:  Comparison  of  first  and  second  order  approximations  of  e'hl  -  1,  which  is  the  form 
of  the  term  responsible  for  phase  and  amplitude  shifts  due  to  interface  height.  The  first  order 
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approximation  is  —  1  fa  i/17,  and  the  second  order  approximation  is  e'h~’  -  1  ih~t  -  |y (h-,  )2. 
(a)  Modulus  of  the  approximation,  (b)  Phase  of  the  approximation  in  radians. 

Figure  16:  RMS  amplitude  error  in  the  comparison  of  perturbation  and  finite  difference  derived 
seismograms.  RMS  amplitude  error  is  measured  after  each  perturbation  method  seismogram 
has  been  shifted  in  time  relative  to  the  complementary  finite  difference  seismogram  such  that 
the  RMS  amplitude  error  is  minimized.  The  error  for  each  model  is  the  average  of  the  RMS 
errors  of  the  seismograms  in  the  model.  This  mean  error  is  plotted  against  RMS  height  for  two 
values  of  RMS  slope,  and  against  slope  for  a  fixed  RMS  height.  Interface  height  is  in  units  of 
Si  wavelengths.  The  data  points  are  labeled  with  their  model  names. 

Figure  17:  Three-dimensional  rough  interface  scattering  model.  The  auto-correlation  function  of  the 
interface  is  a  two-dimensional  Gaussian  with  x  and  y  correlation  lengths  of  2.4  Si  wavelengths 
(0.19  km),  an  RMS  height  of  0.125  Si  wavelengths  (0.010  km),  and  an  RMS  slope  of  0.30. 
Contours  and  axes  are  labeled  in  kilometers. 

Figure  18:  Transmission  coefficients  for  the  model  in  Figure  17.  The  source  is  an  18  Hz  SY  plane 
wave  at  normal  incidence  with  particle  motion  in  the  x-direction. 

figure  19:  Cross-sections  of  the  reflection  and  transmission  coefficients  in  Figure  18  along  10°.  45°. 
and  80°azimuths.  PU,  SU,  and  TU  and  PD,  SD,  and  TD  are  the  up-  and  down-going  P,  SV. 
and  SH  scattering  coefficients. 

Figure  20:  Transmission  scattering  kernels  for  the  model  in  Figure  17. 

Figure  21:  Transmission  coefficients  for  the  model  in  Figure  17.  The  source  is  an  18  Hz  P  plane 
wave  at  normal  incidence. 

I  igure  22:  Cross-sections  of  the  reflection  and  transmission  coefficients  in  Figure  21  along  10°.  45°, 
and  80°azimuths.  PU,  SU,  and  TU  and  PD,  SD,  and  TD  are  the  up-  and  down-going  P,  SY. 
and  SH  scattering  coefficients. 
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ABSTRACT 


A  theoretical  and  computational  tool  is  established  for  the  numerical  evaluation  of  spectral  am¬ 
plitudes  of  point  sources  in  multilayered  anisotropic  crustal  models.  First,  we  obtain  an  explicit 
representation  of  the  spectral  elastodynamic  Green’s  tensor  in  general  anisotropic  media  as  a  sum 
of  three  integrals  over  the  corresponding  three  slowness  surfaces.  The  multidimensional  stationary 
phase  principle  is  then  applied  to  derive  an  asymptotic  approximation  at  the  fa r- Field .  The  theory 
is  applied  to  azimuthally-isotropic  media  for  which  we  offer  an  alternative  representation  of  the 
Green's  tensor  and  its  ensuing  displacements  fields  in  the  form  of  an  exact  Hankel  transform  over 
the  horizontal  wave-number  variable.  The  total  field  is  specified  here  in  terms  of  two  potentials:  an 
SH  potential  and  a  mixed  quasi  transverse-quasi  longitudinal  potential,  both  of  which  assume  the 
role  of  two  scalar  Green’s  functions.  The  availability  of  the  Green’s  tensor  in  analytical  form  enables 
one  to  obtain  readily  numerical  solutions  for  a  wide  selection  of  media  and  sources.  Far-field  spec¬ 
tral  amplitude  radiation  patterns  for  various  sources  and  displacement  components  are  calculated 
and  displayed.  It  is  shown  that  the  radiation  field  of  an  explosion  has  the  following  new  features: 
(1)  Quasi-transverse  waves  are  created  with  four  and  eight  lobe  patterns;  (2)  Quasi-longitudinal 
waves  are  created  for  the  collatitudinal  displacement  with  four  lobe  patterns:  (3)  The  energy  ratio 
SV/ P  may  reach  the  value  of  20  for  more  than  50%  of  the  azimuths  in  crustal  structures  such  as  tuff 
and  shales;  (4)  Radiation  patterns  for  vertical  shear  waves  are  created  which  are  indistinguishable 
from  corresponding  waves  produced  by  earthquake  faults. 

Haskell-type  matrix  algorithm  is  established  for  a  multilayered  azimuthally-isotropic  half-space 
which  enables  one  to  calculate  body  waves  and  surface  waves  in  real-earth  crustal  models. 


1.  INTRODUCTION 


Measurements  of  seismic  velocity  anisotropy  in  exploration  geophysics  from  travel-times  of  P,  SV 
and  S1I  waves  (Weatherby  et  ai,  1934;  McCollum  and  Snell,  1944;  Ricker,  1953;  White  and  Sen- 
gbush,  1953;  Cholet  and  Richards,  1954;  Uhrig  and  Van  Melle,  1955;  de  Segonzac  and  Laherrere, 
1959;  Richards,  1960;  Gretener,  1961;  Van  der  Stoep,  1966;  Robertson  and  Corrigan,  1983;  White 
et  ai,  1983;  Helbig,  1984;  Banik,  1984;  Leary  and  Henyey,  1985;  Stephen,  1985;  Winterstein,  1986; 
Kerner  et  ai,  1989)  have  disclosed  that  many  rocks  in  sedimentary  basins  exhibit  significant  degrees 
of  anisotropy.  Similar  results  were  obtained  for  the  earth’s  crust  and  deep  mantle  via  earthquake 
seismology  (Gupta,  1973;  Crampin,  1985;  Vinnik  et  al.,  1986;  Kafka  and  Reiter,  1987),  and  labo¬ 
ratory  ultrasonic  velocity  measurements  (Giesel,  1963;  Schock  et  ai,  1974;  Levin,  1979;  Jones  and 
Wang,  1981;  Carlson  et  ai,  1984;  Lo  et  ai,  1986;  Thomsen,  1986;  Rai  and  Hanson,  1988).  The 
conclusions  for  exploration  seismology  (Corrigan,  1989)  are: 

•  Transverse  isotropy1  is  ubiquitous,  being  observed  over  a  wide  range  of  depths,  at  many  scales 
of  length  and  in  many  basins  of  economic  interest. 

•  Azimuthal  anisotropy  is  present  in  some  areas  and  can  significantly  impact  the  quality  of 
shear  wave  data.  Its  utility  as  a  tool  for  detecting  fractural  zones  is  still  under  active  study. 

•  P-wave  anisotropy  is  generally  of  the  order  0-10  percent.  Its  effects  may  often  be  masked  by 
heterogeneity  and  noise. 

'In  exploration  and  earthquake  seismology,  the  name  transverse  isotropy  is  usually  associated  with  a  vertical  axis 
of  symmetry  (horizontal  isotropy  plane)  while  azimuthal  anisotropy  is  associated  with  a  vertical  plane  of  symmetry 
(horizontal  axis  of  symmetry). 
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•  S-wave  anisotropy  is  usually  larger  than  P-wave  anisotropj — it  may  easily  reach  10-30  per¬ 
cent.  Anisotropy  must  be  taken  into  account  for  many  applications  of  S-wave  seismology. 

It  has  also  been  previously  known  that  large  volumes  of  inherently  isotropic  rock  may  be  micro- 
fractured  by  regional  tectonic  stresses,  resulting  in  a  system  of  oriented  fractures  along  principal 
axes  which  in  turn  could  lead  to  global  anisotropy  (Bamford  and  Nunn,  1979;  Leary  et  ai,  1987). 
Mechanical  anisotropy  can  also  result  from  the  presence  of  crystals  of  particular  symmetries  or  pe¬ 
riodic  thin  lamination.  All  these  anisotropy  types  are  susceptible  of  exhibiting  shear  wave  splitting 
(birefringence)  where  the  faster  wave  is  polarized  parallel  to  the  fracture  system,  while  the  slower 
is  polarized  perpendicular  to  the  fracture  system. 

The  existing  theoretical  models  lag  far  behind  the  observations.  We  can  group  it  in  three 
categories: 

1.  Propagation  of  plane- waves  in  homogeneous  anisotropic  structures  (Musgrave,  1954,  1970; 
Fedorov,  1968). 

2.  Radiation  from  point-forces  with  no  explicit  analytical  results  (Yeatts,  1984;  Hanyga,  1984). 

3.  Numerical  evaluation  of  WKBJ  seismograms  (Garmany,  1988;  Singh  and  Chapman,  1988). 

Buchwald  (1959)  derived  an  integral  representation  of  the  field  in  terms  of  triple  Fourier  integrals. 
These  were  estimated  asymptotically  using  the  stationary  phase  approximation,  following  a  method 
of  Ligh thill  (1960).  Yeatts  (1984)  solved  the  same  problem  using  the  Radon  transform.  His  results 
are  presented  in  a  form  that  is  suited  exclusively  for  asymptotic  ray  theory.  It  is  also  suspected 
that  his  final  formal  solution  in  the  frequency  domain  lacks  that  part  which  arises  from  integration 
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over  imaginary  plane-wave  angles  [see  his  equation  (27)]. 

The  present  paper  extends  seismic  source  theory  to  general  anisotropic  media.  It  puts  emphasis 
upon  explicit  exact  and  approximate  analytic  evaluation  of  the  Green’s  tensor.  The  elastodynamic 
fields  of  various  point-sources  are  then  systematically  constructed  from  this  tensor.  In  this  sense,  our 
results  go  beyond  the  works  mentioned  above,  since  we  equipped  seismologists  with  a  theoretical 
tool  by  which  they  can  interpret  observations  of  amplitudes  of  seismic  waves  as  well  as  travel- 
times.  The  theory  was  applied  to  media  with  azimuthal-isotropy:  Here,  due  to  the  complexity 
of  the  analysis,  it  was  necessary  to  provide  two  parallel  theoretical  approaches  through  which  the 
numerical  results  could  be  checked  and  compared.  To  this  end  we  have  presented  the  displacement 
far- field,  once  in  a  form  of  a  stationary-phase  approximation  and  then  alternatively  in  the  form  of 
an  exact  Hankel-Transform  over  the  horizontal  wave-number  variable.  The  results  of  this  paper 
are  ‘computer- ready’  in  the  sense  that  the  exact  expressions  for  the  displacements  due  to  a  variety 
of  dipolar  sources  in  an  homogeneous  layer  are  given.  When  applied  to  a  multilayered  media  with 
proper  boundary  conditions,  Haskell  or  matrizant  type  algorithms,  would  render  an  immediate 
generalization  of  our  results. 
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2.  INTEGRAL  REPRESENTATION  OF  THE  DISPLACEMENT  FIELD 


We  consider  the  equation  of  motion  for  an  elastic  anisotropic  medium  (Ben-Menahem  and  Singh, 
1981) 

rflu-  fftni 

(2.1) 


d2u.i  d2ui 

C'Jkt  p-  Q_  — 


dt 2  3  dxjdxk 

where  u,  =  u,(r,  t)  is  the  i-th  component  of  the  displacement,  /,  =  /,(f,  t)  is  the  i-th  component  of 
the  applied  body  force  per  unit  volume,  p  is  the  uniform  mass  density,  and  c,:kt  are  the  elements 
of  the  uniform  elastic  coefficient  tensor  which  satisfy  the  symmetry  conditions 


^ijkl  —  Cjtkl  —  C ijlk  —  Ckliji  (2.2) 

so  that  only  21  independent  constants  are  involved.  The  suffixes  can  take  the  values  1,2,  or  3,  and 
the  summation  convention  for  repeated  suffixes  is  assumed. 

Let  the  applied  body  force  be  an  arbitrary  point  force  F,(t)  acting  at  the  origin.  Then 

f,(r,t)  =  F,{t)6(r)  ,  (2.3) 

where  S(r)  is  the  three-dimensional  Dirac  delta  function. 

We  assume  that  the  force  and  the  displacement  can  be  represented  by  the  Fourier  integrals 

1  f°° 

Ft(u)=—  dte-**Fi{t),  (2.4) 

Z7T  J  —  oo 

and 

■  (i)7 .  (2.5) 

In  (2.4)  we  have  used  the  notation  F,  to  describe  two  different  functions,  emphasizing  with  the 
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arguments  the  corresponding  domain.  The  same  applies  to  u,  in  (2.5).  Thus 


/oo 

du>etwiFi{ w)  , 

-OO 


u,(r,t)  =  j°°  duj  J  J  iPke'^-^Uifcu)  ,  (2.7) 

where  u  is  the  angular  frequency  and  i  =  \/-T. 

Taking  the  Fourier  transform  of  (2.1),  and  using  (2.4)  and  (2.5),  we  obtain  the  matrix  equation 


( ciijkkikk  -  pu26ij)uj(k,u>)  =1—1  Fi(u)  . 


Letting 


k  =  kh  , 


where  n  is  a  unit  vector  in  the  direction  of  k,  and 


A  Clijk 

Aij  =  — —niTik  , 

P 


(2.10) 


Equation  (2.8)  takes  the  form 


f«  *'-(£)«’ 


(2.11) 


where 


The  formal  solution  to  (2.11)  is 


MtJ  =  pk2  Av-  -  6. 


U'  ~  (27)  M'J  lp] 


(2.12) 


(2.13) 
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We  may  write  (2.13)  in  an  explicit  way  by  means  of  the  spectral  resolution  (otherwise  known  as 
the  “nonion  form  of  a  dyadic”)  of  the  matrix  A  =  [/It;]  ;  that  is. 

A,J  ---  •£  ClA 17'  .  (2.M) 

m~  1 

where  C ^  is  the  m-th  eigenvalue  of  A  and  id)”1'  is  the  m-th  projection  matrix.  Since  A  is  real, 
symmetric,  and  positive  definite,  its  eigenvalues  are  real  and  positive. 

The  eigenvalue  equation  for  A  is  the  Christoffel  equation  of  plane-wave  theory 

AX}T/m)  =  -  (2.15) 

wherein  Cm  is  the  phase  velocity  of  type-m  waves  and  is  the  corresponding  displacement 

polarization  vector  (Fedorov,  1968). 

The  eigenvectors  can  be  taken  as  orthonormal, 

r,(p)T,(?)  =  S„  ,  (2.16) 


and  the  projection  matrices,  for  symmetric  A ,  can  be  expressed  as 


AXJ(m)  =  . 


Consequently,  the  eigenvalues  ( (3m )  of  the  matrix  M_  are 


Pm  =  pk 2 


Cl 


-  (f); 


,  m  =  1,2,3 


and  the  eigenvectors  of  M.  are  the  same  as  those  of  A  . 
Therefore,  we  write  for  the  inverse  of  the  matrix  M 


Mi,~'  =  t  W-Avml- 

m  =  1 


(2.17) 


(2.18) 


(2.19) 
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The  solution  in  the  (A;,u;)-domain,  (2.13),  is  then  written  explicitly  as 


where  we  have  emphasized  the  dependence  of  and  of  Cm  on  h  . 


(2.20) 
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3.  THE  SPECTRAL  GREEN’S  TENSOR  AND  ITS  SLOWNESS-SURFACE 

INTEGRAL 

Equation  (2.20)  can  be  transformed  into  the  (f,  u)-domain  by  standard  procedures.  Skipping  the 
details  of  this  lengthy  operation,  the  end  result  for  the  outgoing  field  is 

u,(f»  =  E  /  dSlAt]WC-3H(r.  h)e~Mlc”'  ,  (3.1) 

'  m=  1  ^ 

where  dfl  is  the  solid  angle  element  and  H  is  the  Heaviside  unit-step  function.  We  write 

Ui(r,u)  =  Gij(r,u)Fj(u),  (3.2) 

where  the  Green’s  tensor  Gij(r,w)  is  given  explicitly  by 

E  J  dSlAi/m)C-3H(r.  h)e~M'c -  .  (3.3) 

'  m=  1 

If  fl  is  chosen  as  the  unit  sphere,  the  representation  (3.3)  can  be  written  in  the  equivalent 
expression 

=  Cd*  L  d^^A^m)(h)C'3(h)H(f-  n)  • 

'  m=l 

‘  exP  t  77-7-rT  [sin  cos  0  +  y  sin  ^)  +  zcosO]  l  ,  (3.4) 

lCm(n)  J 

where  C±  in  (3.4)  is  the  Weyl  contour  {0  <  6  <  x}  U  {0  =  tt/2  ±  it,  (t  >  0)}  specified  in  Figure 
1. 

In  the  particular  case  of  an  isotropic  medium,  the  eigenvalues  C^h)  and  the  projection  oper¬ 
ators  A,/m)(h)  assume  the  form 

C?  =  cl  =  C|  =  t  3  (3.5) 
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we  obtain  from  (3.9) 


(3-11) 


which  is  the  well-known  expression  for  the  Green’s  tensor  in  an  isotropic  medium  (Ben-Menahem 
and  Singh,  1981). 


63 


The  Slowness  Surface 


The  surface  whose  position  vector  has  the  magnitude  1  /Cm  is  usually  called  the  slowness  surface 
(Synge,  1957). 

Let  s  be  the  position  vector  of  a  point  on  the  slowness  surface  Tm  of  type-m  waves,  then 


n 


s  = 


(3.12) 


A  point  P  on  the  slowness  surface  with  coordinates  s,  represents  a  wave  travelling  in  the 
direction  OP  with  velocity  l/|s|  .  The  slowness  surface  is  given  by  Cm  =  1,  m  =  1,2,  or  3,  and 
for  each  m  ,  Cm(s)  represents  a  real  sheet  Tm.  Since  Cm(s)  >  0  depends  continuously  on  h,  where 
s  6  Tm  if  s  =  h/Cm,  it  follows  that  the  rm  are  three  closed  surfaces  enclosing  the  origin. 

Since  the  normal  to  the  slowness  surface  coincides  with  the  direction  of  the  energy  flux,  we  can 
introduce  the  unit  normal  vector  N  (Musgrave,  1970), 


W  m  , 


(3.13) 


where  tlnm)  is  the  group  velocity  vector 


cijkl  rji(m) 

pCm  3 


rp(m) 


Til  = 


—  c'ikl  a(m) 


-Si, 


(3.14) 


and  wm  is  its  magnitude. 

Next,  we  express  the  Green’s  tensor  in  (3.3)  as  an  integration  over  the  slowness  surface.  Using 
(3.12)  we  modify  (3.3)  as  follows.  First,  from  (2.10),  we  have  on  Tm 


A,J(s)  =  C~2A,J(n). 


(3.15) 
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By  comparing  (3.15)  with  the  eigenvalue  equation  (2.15),  we  find  that  on  Tm  the  eigenvalues 
are  equal  to  unity. 

The  second  change  involves  the  surface  element  on  Tm 

surface  element  =  sidai  ,  (3.16) 


where  doi  is  a  component  of  the  outwardly  directed  surface  element  of  Tm  and  s/  is  as  in  (3.12). 
On  the  slowness  surface 

dfft  =  Nida  .  (3.17) 


Then,  using  (3.12)  and  (3.13),  we  have 


•  (m)  d(T 

sidai  =  niw)  - 


da 


(3.18) 


where  we  have  used  the  fact  that  (Lighthill,  1960), 


(m) 

niw)  =  Cm  , 


(3.19) 


i.e.,  the  velocity  of  energy  propagation  normal  to  the  wave  fronts  is  equal  to  the  phase  velocity. 
Using  the  above  results  in  (3.3)  we  obtain  for  the  Green’s  tensor 

=8 S  ^  /  d°  s)e~'“?*  ,  (3.20) 

”  m=l  "'O'1 

where  the  integration  is  now  over  the  slowness  surface  Tm  of  type-m  waves. 

In  this  way,  instead  of  specifying  a  wave  normal  h  and  a  velocity  Cm  of  the  wave  parallel  to  n, 
we  employ  the  slowness  vector  s  =  h/Cm. 

Equation  (3.20)  for  the  Green’s  tensor  can  be  written  as 


Gt](r,u)  =  Gj(f,W)  +  G^(r,w)  +  Gg"(f,w)  • 


(3.21) 


m=  1 


m  —  2 


m= 3 
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This  separation  of  the  Green's  tensor  in  three  terms  is  merely  mathematical  and  is  not  necessarily 
associated  with  true  P,  SV,  or  SH  motions.  However,  for  instance,  in  homogeneous  media  with 
hexagonal  (azimuthal  isotropy)  and  cubic  symmetries,  it  is  easy  to  show  that  true  SH  motion 
(understood  as  divu  =  0)  is  possible,  and  consequently,  the  concept  of  SH  waves  can  still  be 
maintained. 

In  this  case,  the  SH  part  of  the  total  Green’s  tensor  is 

GjV.w)  =  /r  •  (3.22) 
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4.  AZIMUTHALLY  ISOTROPIC  MEDIA 


4.1  SH  Motion 


In  this  section  we  apply  the  previously  developed  theory  to  the  case  of  the  azimuthally  isotropic 
medium,  where  we  only  have  five  independent  elastic  constants 


Cnn  —  C2222  =  cn 

C3333  =  C33 

C1122  =  C12 

(4.1) 

Cl  133  =  C2233  =  C13 

C1313  =  C2323  =  C44 

C1212  =  C66  =  (Cn  -  C12)/ 2  . 

The  remaining  elastic  constants  are  equal  to  zero  if  not  related  to  the  above  by  the  basic  symmetry 
relations  (2.2). 

In  an  azimuthally  isotropic  medium  we  have  for  the  displacement  polarization  vector  of  SH 
waves 

I-,) 


f(3)(ii  =  (s\  +  sir1'2 


Si 


(4.2) 


\ 


0 


) 


Therefore,  the  corresponding  projection  operator  takes  the  form 


( 


4(3)(^  =  (5?  +  sir1 


-sis2 

0 


-sis^  0 
sf  0 
0  0 


(4.3) 


67 


For  these  purely  transverse  waves,  the  phase  velocity  (eigen\alue)  can  be  calculated  from  (2.15) 


fa(.s-)  = 


<Y>« 


P  +  (c66  —  C44  )«? 


1/2 


(l-l) 


The  group  velocity  vector  for  SII  waves  can  be  obtained  from  (3.14), 


=  —  (siSi  +  s2h)  +  —  53«3 


c66 


C44 


We  proceed  to  obtain  the  equation  of  the  slowness  surface.  From  (3.12)  we  can  write 


|s|2  =  si  +  s]  +  sl  =  C32  . 


(4.5) 


(4.6) 


Substituting  (4.4)  into  (4.6)  we  find 


— w  +  »p  +  —  4  =  1 

P  P 


(4.7) 


which  represents  the  equation  of  the  slowness  surface  for  SH  waves.  This  surface  is  an  ellipsoid, 
and  in  the  case  of  isotropy,  (c44  =  c66  =  ;r),  reduces  to  a  sphere  of  radius  \fpjji. 

To  perform  the  integration  in  (3.22)  over  the  slowness  surface,  we  change  to  prolate  spheroidal 
coordinates  (77,  9,0),  where 

x  =  a  s i n h 77  sin  9  cos  0 

y  =  o  sinlwy  sin#  sind>  (4-s) 

z  =  a  cosh  7  cos  9 

with  0  <  7  <  oe, 0  <  9  <  7T,  and  0  <  <  27r.  The  corresponding  met  ic  coefficients  are 


Gw  =  Ooe  =  «2(sinh27  +  sin 2  9) 
g.;,#  =  a2  sinh2 7  sin2  9 


[  1.0) 
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For  rj  =  constant  we  have  the  prolate  spheroids 


a2sinh2r/(2:2  +  y2)+ a^cosh2/2  "  1 


(4.10) 


Thus,  with  the  use  of  prolate  spheroidal  coordinates,  the  integration  over  the  slowness  surface  is 
carried  out  at  rj  =  constant.  ^.Frorn  (4.7)  and  (4.10)  we  find  that 


■  £44  V  66  /  . 


(4.11) 


t)  =  tgh  1  f  — )  ,  (4.12) 

\c66/ 

which  is  the  equation  of  the  slowness  surface  in  prolate  spheroidal  coordinates. 

For  the  area  element  on  the  slowness  surface  we  have 

do  =  (geeg<t><t>)Xl2d9d<p  =  —  sin0  fcos2  9  +  —  sin2  §)  d9d<i>  ,  (4.13) 

C66  V  C44  J 

whereas  the  magnitude  of  the  group  velocity  on  the  slowness  surface  is 


-  (t)  (■ 


1/2  /  _  \  1/2 
(cos2<Si-f  —  sin2  0  J 


Also,  the  projection  matrix  takes  the  form 


sin2d>  -sind>cos</>  0 
\4>)  —  -sind>cos</>  cos2  A  0 


Therefore,  from  (3.22),  (4.13),  and  (4.1  4)  we  have  for  the  Sli  part  of  the  Green's  tonsi 


5v.*) = )1/2  r  d0AtJ^(O)  1 

»nZCr.c  \C44J  J0  Jr± 


(If)  sili  9  //(  r  ■  s)t  ,~r 


Equation  (4.16)  represents  essentially  an  integration  over  the  complex  unit  sphere  as  we  had  in 


(3.4).  Since  the  direction  of  the  polar  axis  is  arbitrary  and  fc±  -  f0/,2±'~>  we  can  omit  the  step 


function  Il(r-s).  Consequently, 


C  d>A^L a 


(4.17) 


The  projection  operator  in  (4.15)  can  be  rewritten  as 


cos2  d>  sin<i>cosd>  0 
—  L  ~  sin<pcos0  sin2<p  0 


0  1 


where  I  is  the  unit  matrix. 


We  recast  the  exponential  in  (4.17)  in  the  form 
(~'~r  5  =  exp  |  —  ^  (x  cos<p  +  ys\n(p)s\nd  +  ^  —  ^  xcos# 


(4.19) 


where 


-v  /  /  \  /  /  ^*66 

r  =  (x,y,x  ),  x  =  zl  — 

\  Ca 4 


n  —  (cos<psin#,sind>sin0,cos#)  .  (4.21) 

Inserting  (4.19)  into  (4.17),  one  shows  that  (4.18)  can  be  put  in  the  symbolic  operational  form 

4(3)(0)  =7  ~eze2  ~  .  (4.22) 


where 


rAu  =  <±<± '/  x-  I  >  ...  -  '  v  a  i-^y  f  -V/(:.A)r/A  . 


7f) 


and 

+  *  +  +  (426) 

Explicitly, 

Gsh  (r,u)  =  e^giz,  A)  +  ( eAeA  -  e^)-^  J  Ag{z,A)dA  ,  (4.27) 

which  agrees  with  the  result  obtained  by  Ben-Menahem  (1989).  We  also  note  that  the  integral  in 
(4.27)  renders 

/AS(J'A)<,A  "  4 '  ‘4'28) 

In  Appendix  A  we  give  certain  differential  geometry  results  for  the  slowness  surface  based  upon 
a  general  parameterization.  These  results  will  be  useful  for  the  evaluation  of  the  Green’s  tensor  in 
systems  with  lower  symmetries. 


4.2  Quasi-Longitudinal  and  Quasi-Transverse  Motions 

We  represent  the  explicit  Green’s  tensor  in  terms  of  a  single  field  potential  for  the  coupled  P-SV 
field.  The  details  of  this  operation  are  sidetracked  to  Appendix  B. 
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We  proceed  to  solve  the  simultaneous  equations  (B.8)  and  (B.9).  To  this  end,  we  follow  the 


notation 

Cll  C33  C13  +  C44  C66  C  44 

ai  =  —  ;  a2  =  —  ,  a 3  -  - ,  a4  —  —  ,  a 5  -  — 

P  P  P  P  P 

and  put  (B.8)  and  (B.9)  in  the  form 

(a5V<  +  a2d]  +  w2)r  +  a3d2®  =  -fe  ■  {ezdzS(f-  f0)}  , 
a3V2 r  +  (ai  V2  +  asd2  +  oj2)Q  =  -fe  ■  Vt6(r-  r0)  . 

Denote 

L\  =  asV2  +  a2d]  +  n>2  , 
i/2  =  aiV2  -f  a$d2z  +  cj2  , 


(4.29) 


(4.30) 

(4.31) 


(4.32) 

(4.33) 


and  assume  a  representation  of  the  auxiliary  potentials  0  and  V  in  terms  of  a  single  fundamental 
field  potential 


Fe-[ezdza3V2-VtLx}V  , 

(4.34) 

Fe-[-ezdzL2  +  a3Vtd2z]V  . 

(4.35) 

A  substitution  of  (4.34)  and  (4.35)  into  (4.30)  and  (4.31)  reveals  that  is  the  scalar  Green’s 
function  of  the  equation. 

(-LiL2  +  a23d2Vf)^  =  —6(f—  f0)  .  (4.36) 

In  the  derivation  of  (4.36),  we  were  aided  by  the  fact  that  the  linear  operations  L\,  L2,  dz,  V2, 
and  V,,  commute.  Since  0  and  f  are  expressible  in  terms  of  ty,  so  can  the  displacements.  Indeed, 
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using  (B.10)  we  have 


up-SV  =  Vj{Vt  20}  +  ez  J  Tdz  =Gp-sv  Fe  .  (4.37) 

We  put  here  the  expressions  for  0,T  from  (4.34)  and  (4.35).  Some  rearrangements  are  necessary 
in  order  to  arrive  at  the  explicit  expression  for  the  P-SV  part  of  the  spectral  Green’s  tensor 

Gp-sv  {F\f0;u)  =  -  a3W  -  ezezL3 - 4tT'^4  ^  .  (4.38) 

P  v<  . 

with 


—  a  1^?  +  (<13  +  Qs)d2  +  , 

(4.39) 

l4 

=  (&3  +  +  to2  , 

(4.40) 

vv 

=  ('9* + *■£)  (** + s-l) 

(4.41) 

=  ezezd2z  +  eAeAdA  +  e<t>e<t>—dA  +  {ezeA  +  eAez)dAz  . 

(4.42) 

Clearly,  0$  -  0  due  to  symmetry  and  differentiation  is  w  ■  r  ■  t  field  coordinates.  Nevertheless,  one 
must  note  that  the  operators  L4  and  V(Vf(Vj"2)  [see  ( B .32)]  do  not  commute,  i.e.,  VtVf( V~2T4)  ^ 
The  total  Green’s  tensor  is  the  sum  of  the  expressions  in  (4.38)  and  (B.36).  Explicitly, 

4t rp  G  (f/f0]u)  =  47 ip  G„a  -47r p  Go 

s  $ 

-  e2ez  - +  47tL3{$}  -  eAeA  — - —  +  47rT4{^}  ,  (4.43) 

IVSHO-S  VSH 

4?rp  Gtso=  — - —  I  5  +  47Ta3VV^,  (4.44) 

VSHd  5 

4? rp  Go=  (e*e*  -  eAeA)-^r  j  A  +  47rZ,4{^}l  dA  .  (4.45) 

A2  J  l  Vsna s  J 
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The  arbitrary  partition  of  the  total  Green’s  tensor  in  Eq.  (4.43)  is  meaningful.  In  the  isotropic 

limit  ( tjsh  —  1>  ai  —  <*2  =  ,  04  =  a5  =  03  =  one  finds  that  Go=  0,  Giso  goes  into  the 

Green’s  tensor  for  unbounded  isotropic  media  and  G=Giso- 

In  (4.43),  G  is  symmetric  in  its  indices  and  also  invariant  to  an  exchange  of  rand  fo-  We  shall 
♦— * 

name  Go  the  singularity  tensor  of  the  medium.  We  also  show  in  the  next  section  that  the  potential 
splits  naturally  into  a  quasi-longitudinal  part  $ p  and  a  quasi-transverse  part  such  that 


¥  =  Vs  ~ 


(4.46) 


Accordingly,  we  effect  a  sub-partitioning  of  the  singularity  tensor  into  three  constituents 

«  _SV  —P 

Go-Go  +  G  0  —  Go  1 


(4.47) 


where 


~SH 

PG  0  =  (e*e*  -  eAeA) 


t  e-^snR 

^irksHVsHas  A2 

PGo=  -  eAeA)^2  J  AT4{4';(A,a)}dA  . 

U  =  p,sv) 

We  shall  find  it  useful  for  a  later  development  of  the  theory  to  recast  (4.43)  as 


(4.48) 

(4.49) 


G=Gsh  +  Gsv  -  Gp 


(4.50) 


where 


PGp  =  -a  xeze2 
+  a\e&e& 


V?  +  —d;  +  k 


f  03 

l  o\ 


d2A*p- 


+  a3(ezeA  +  eAez)d'ii9p 


+  a-ld)  +  kl 

Gl  0\ 


1  ~P 

+  03^6^,—  db'S/j,  4-  p  Go  > 


(4.51) 
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and 


pGsv  =  -o.se2e2 


^v?  +  d2z  +  kl 

«S 


Vs  +  a3(ezeA  +  e&ex)d\zV  s 


+  ase  aCaI— 
I  ^5 


«3  +  05  2  .  a2 


v?  +  ^ d\  +  k 2 

05 


05 


.sv 


+  a.3e<pe<p—d£,'Sfs  —  pG  o 


It  can  be  shown  that  the  singularity  tensor  has  the  following  properties: 


1.  In  the  absolute  far-field  (A  — ►  oo,  z  — ►  oo) 


Go=  O  (  £7  )  . 


2.  In  the  axial  far-field  (A  =  0,  z  -»  oo) 

«S//  /  1  \  /1\ 

Go  =  O  (j^J  (exactly)  ;  G0=  O  ^-J  (non-singular),  j  =  P ,  SV 


3.  In  the  isotropic  limit 


Go 


(~P  «SH  ~SK\ 

=  0  [Go  —  0,  Go  =  —  Go  J 


Developing  the  operators  in  (4.51)  and  (4.52),  .ve  find  a  simplified  form  of  the  P-SV  Green’s 

(4.38) 

Gp-sv  (r|rb; w)  =  -  {-eze*[aiV?  +  ahd]  +  u2}$  +  a3(e2eA  +  eAe2)d2AzV 
P  1 

-  +  [a2d]  +  w2)^  J  A'&dA 

-  eAeA  a5dl*  +  (a2d2  +  u;2)  (v  -  J  AtPdA^  |  . 

Thus,  in  the  absolute  far-field,  we  have  the  approximation 


P  Gp~  -o ie,e; 


(d2A  +  +  k lj  yp  +  a3(e2eA  +  cAc2)dAzVf 


(4.52) 


(4.53) 


(4.54) 


(4.55) 

tensor 


(4.56) 
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-  ( a±dl  +  +  kl)  *„  (4.57) 

Vai  ax  / 

p  Gsv«  -a5ezez  — d |  +  d2  +  fc*l  +  a3(e2eA  +  eAe2)^z*i 
1^5  J 

-  (»i  +  |a,2  +  it;)  *.  •  (4.58) 

iFrom  (4.55)  and  (4.58)  we  may  derive  in  a  straightforward  manner,  the  P-SV  far- field,  of  point 
forces  in  terms  of  the  yet  undetermined  potentials  'fp  and 
I.  Force  in  the  z-direction  (axis  of  symmetry) 


„  F(u>)  «  . 

u  — - G  ‘&zi 

P 


F  p  p  \  a  i  / 

p  p 

»„  =  +  *?  +  «)*.. 
p  p  \05  / 

^(W)  ~  -  -  ^(<*>)  »2  iTl 

=  - -  Gsv ■  cz€a  = - a3aAiw,. 

P  P 

II.  Force  in  the  x-direction  (plane  of  isotropy;  ex  •  e&  =  co6  <f>) 


(4.59) 


..  F(u)  ~  . 

U  =  — ^ -  G  -Cxr 

P 

jFYw)  ~  .  .  F(u>)  ,fl2  lTf 

Upz  = - -  Gp'-  exez  =  —  - - ~a3  cos^^p, 

P  P 

F(<jj)  «  .  -  F(w)  _  JL  r®5  a2  ,  °2  fl2  ,  1.2I  ,T, 

UpA  = - L  Gp •  creA  =  - °i  cos  <t>  ~°A  +  ~d*  +  kr  *P’ 

p  p  Lai  °i 

u«  =  Gsv:  exe*  =  ^^a3cos 4>d\zV t, 

P  P 

«,A  =  Gsv:  ereA  =  a5cos<f>  d\  +  —  d]  +  k]  <Sl,  .  (4.60) 

p  p  as 

Expressions  (4.57)-(4.60)  can  also  be  used  in  the  axial  far-field  with  insignificant  error  for  most 


crustal  models. 


4.3  The  Fundamental  Field  Potential 


In  the  previous  section,  the  non-SH  field  was  found  to  depend  on  a  single  potential.  Consequently, 
we  have  given  the  curvilinear  components  of  the  Green’s  tensor  in  terms  of  this  potential.  It 
remains,  therefore,  to  solve  (4.36).  This  is  accomplished  by  means  of  the  three-dimensional  spatial 


Fourier  integral.  Assume  the  feasibility  of  the  representation 


*  =  £///_”  7  •*•**«>  ***,. 


(4.61) 


and  the  known  identity 


6(f)  =  _L  J  J  J°°  e-<°*+0v+-r*)  dadpdy  . 


(4.62) 


Inserting  these  in  Eq.  (4.36)  it  is  found  that 


Y(a,/3,7;u;)  =  K72  +  m(a2  +  /l2)  -  w2][a5(a2  +  (?)  +  a272  -  a;2]  -  a272(a2  +  01)  ■  (4.63) 


We  assume  here  without  loss  of  generality  that  r0  =  0,  since  it  can  be  reinstalled  in  the  final  result. 
We  make  a  change  of  variable  from  (a,  /3,7)  to  spherical  coordinates  in  wave-number  space 

a  =  —  sin  0  cos  ^  ;  /?  =  —  sin  0  sin  0  ;  7  =  —  cos  0, 

000 


a2  +  p2  +  72  =  ^r- 


(4.64) 


The  Jacobian  of  this  transformation  is  J  =  j^-sin^|.  Some  algebra  is  required  to  show  that  Y  in 
(4.63)  is  indeed  factorable 


Y  =  £4  (C2  -  C2)(C2  -  Cly)  , 


(4.65) 


where 


2  m  +  y/M  2  m  ~ 

Cp  ~ - 2 -  ’  °sv  - - 2 - 
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m  =  a5  +  d2  cos2  9  -f  ax  sin2  9  , 


(4.66) 


M  =  [(ai  -  a5)sin2  5  + 


l 


2Q3 


ai  —  as 

2  4)  ^ y. .  j*.  \  /3l2  _j_  2  « 


J  cos2  9 

2 

~  4a| 

/ 

a2  -  (ai  -  a5)(a2  ~  as) 
(d  -  a5)2 


cos 
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=  ((d!  -  d5)  sin2  9  -  (a2  -  a5)  cos2  5]2  +  03  sin2  25  . 


Note  that  Cp  and  Csv  are  the  plane-wave  phase-velocities  of  the  respective  quasi-longitudinal  and 
quasi-transverse  waves  (e.g.,  Fedorov,  1968;  Musgrave,  1970).  Only  in  the  isotropic  limit  do  these 
velocities  signify  a  true  longitudinal  and  transverse  motions,  respectively. 

It  is  apparent  from  (4.66)  that  for  media  in  which 


a\  =  (di  -  a5)(a2  -  as)  , 


(4.67) 


the  phase  velocities  assume  the  simpler  form 


C2  =  (ax  sin2  9  +  a2  cos2  5)  ;  C|v  =  a,  . 


(4.68) 


Such  media  are  known  as  ellipsoidally-anisotropic  (Gassinan,  1964).  The  P-wave  fronts  are  ellipsoids 
of  revolution.  Thomsen  (1986)  showed  that  it  is  a  special  case  of  weak  anisotropy,  obtained  by  a 
Taylor-series  expansion  of  the  phase-velocity  functions  in  (4.66)  and  (B.21).  Keeping  the  first  few 
terms  therein,  one  obtains: 


where 


Cp 

Csv 

Csh 


•v/aJt  1  +  6 sin2  9 cos2  9  +  esin4  9]  +  0(e2;£2), 


y/a$ 

1  +  — (e  -  5)sin2  0co62  9 

+  0 

r(e " 6) 

L  a5 

J 

r 

La5 

1  +  g  (Vsh  ~  1)  s*n2  ^ 

+  0 

(2k^) 

ai  ~  fl2  .  ,  a2  ~  (°2  ~  as)2  .  02,  _  t\  _  I  ( _L  _  1)  ai  ~  °5 

2a2  ’  ~  2a2(a2  —  a5)  a5'f  2  \  y?  J  a2  -  a5  ' 


(4.69) 

(4.70) 

(4.71) 


(4.72) 
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The  case  e  =  6  corresponds  to  ellipsoidal  anisotropy.  It  is  appropriate  at  this  point  to  introduce 
some  measure  of  anisotropy  in  the  form  of  percentage  deviation  from  isotropy.  We  arbitrarily  define 
this  deviation  in  terms  of  the  parameter  e. 


Inserting  (4.64)  and  (4.65)  into  (4.61),  the  fundamental  integral  becomes 


*  =  _JL_  r 

8 tt3w  Jo 


_ _  _ [  d9 sin 6e~'&ZCOBS  F*  e-‘rsin®(xcos^+»'s,n'*)  di 

[C2  -  ctmc*  -  Clv(0)]  Jc±  Jo  * 


(4.73) 


Since  7  =  J (q2  +  /32)  -  ^7,  where  (a2  +  /32)  vary  from  zero  to  infinity,  7  can  obtain  pure  imaginary 
values.  The  variation  of  7  is  then  mapped  onto  a  corresponding  variation  of  0  along  either  C+  or 
C~  (Fig.  1),  depending  on  the  sign  of  z,  to  secure  finite  values  of  e~nz.  For  all  values  of  6,  C2(0) 
and  C|v(0)  are  real.  In  (4.73),  is  the  modified  Weyl  contour,  as  explained  in  Appendix  B. 
The  integral  over  <j>  is  immediate  and  equals  2irJo  (fA  sin  0),  where  J0  is  the  Bessel  function  of  the 
first  kind  and  order  zero.  Then,  the  integral  over  C  is  evaluated  through  the  method  of  residues  at 


the  contributing  poles  C  =  Cp  and  C  =  Csv ■  Using 

J- cc  (C*  -  C$  +  w)(C2  -  C2SV  +  10)  C2-C2SV  lCsvC 


->(c^Asinf0 


1  —ijSj-zcos  Q 

— -  e  cp 


(4.74) 


The  net  result  of  the  two  integrations  is 

lTf  _  f  ddsind  /  (^Asinfl)  -t^-zcoiS  ^(^Asin^)  _ t»_ 

4xu  Jc±  C2(0)-C2sv(0)  |  CSv{0)  e  SV  Cp(9)  6 


(4.75) 


In  the  isotropic  limit  Cv  — *  a,  Csv  -*  b  and  the  classical  Sommerfeld-Weyl  integral  (a  =  a  or  b) 


— 7 — —  =  /  d0sin  9Jo  ( — A  sin  fl')  e  'fzcos^ 

-ikaR  Jc±  V  <t  J 


(4.76) 
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secures  the  reduction  of  ®  to  the  isotropic  limit.  Note  that  the  first  integral  on  the  r.h.s.  of 
(4.75)  represents  the  quasi-transverse  part  of  the  field  while  the  second  integral  is  due  to  the  quasi¬ 
longitudinal  motion.  The  explicit  evaluation  of  $ ,  and  Vp  via  (4.75)  can  be  facilitated  by  reducing 
it  to  a  single  integral.  Indeed,  starting  from  (4.61)  and  changing  variables  and  coordinates  to 
(Fig.  2). 

a  =  k  cos<f>,  f3  =  k  sin<£,  a;  =  Acos<£,  j/  =  Asin<£;  7  =  s,  (4.77) 

one  arrives  at  the  representation 

1  f-oo  foo  r2ir  _  , 

*  =  ^5  J  dse~'aZ  J  kdk  J  d4>e-'k*coal+-*'>  Y~\k\s2) 

1  r°°  ro o 

=  4^  Jo  J°(kA')kdk  J  dae~tsz  Y~1(k2  ,s2)ds  ,  (4.78) 

where 

Y  =  d^s4  +  d2s2  +  do, 


d\  =  02a5, 

d2  =  a2ax{k2  -  k2v)  +  a\{k2  -  k2)  -  a\k2  ,  (4.79) 

do  =  a5ai(fc2  -  kl)(k2  -  k2)  . 

The  integral  over  s  is  evaluated  by  taking  the  residues  at  its  s-poles.  The  final  result  has  the  form 
of  a  Hankel  transform 


»(A,*) 


1  J0(kA)kdk  e-UISi(fc) 

87ra2o5  Jo  VAk*  +  Bk 2  +  E 2  [  S2(/fc)  5i(Ar) 

Vs  -  Vp, 


(4.80) 
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5?  =  A;2  -  ^rj  =  -£**((1  +  ^)+  nA:2  +  v'aFJbFTS7, 

Sl  =  k2-c^Tk)=  ~\kl  (x  +  ftt)  +  n*2  ’  vMA:«  +  Bk 1  +  £*> 


A  =  n2  -  m, 


(4.81) 


v  / 

B  =  (m  —  n)  (l  +  fJJ-)  A:2, 

£2  =  [fl  +  H2  -  r4~l  -  £2  "  4<4£  >  0 

4  P  I  \  C44  y  C33C44 J  ’  — 

m  -  £11.  •  n  _£u _ £13-  £li  _  £11 

m  “  C33  ’  2C44  2c33  C„  C33  ‘ 

The  representation  of  the  Green’s  tensor  as  a  Hankel  transform  can  be  straightforwardly  obtained 
by  substituting  from  (4.80)  into  (4.56). 

In  the  isotropic  limit 

m  =  n  =  1  ,  A  =  B  =  0,  E  =  i(A:2  -  A:2)  ,  St  =  y/k*  -  fc2  ,  S2  =  \jk2  -  k2  ,  (4.82) 

and  the  Sommerfeld  integral 


c—ikaR  roo 


roo 

=  /  J0(kA)e 

JO 


-|*lv**-*f . 


(4.83) 


guarantees  that  (4.80)  falls  back  on  the  isotropic  limit.  Otherwise,  the  integral  in  (4.80)  cannot 
be  evaluated  in  closed  analytical  form  except  for  special  cases.  The  integral  representation  of  the 
potential  in  (4.80)  is  valid  for  all  values  of  the  elastic  parameters  and  its  numerical  evaluation  is 
the  subject  of  a  sequel  to  the  present  paper.  Another  advantage  of  (4.80)  is  its  suitability  for  the 
numerical  calculation  of  displacement  fields  of  sources  in  multilayered  azimuthally  isotropic  media, 
using  the  known  Haskell-Thomson  matrix  algorithm. 
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Fields  of  Point  Forces 


I.  Localized  force  in  the  direction  of  the  symmetry  axis 

This  force  does  not  excite  SH  waves.  The  displacements  are  given  by  u  =  F(u)  G  e2.  Starting 
from  (4.57)  and  (4.58),  we  obtain  for  a  force  in  the  ^-direction 

-sSb 

The  error  term  in  all  expressions  is  O  •  Radial  and  collatitudinal  motions  are  obtained 

via  the  relations 


[a3Vd2  -  e2(w2  +  aiV2  +  (a3  +  a5)d2)  ^  .  (4.84) 


urp  =  uzpcos9  +  u&pSinO  ;  ugp  =  -uzpsinff  -f  u&pcos8, 


(4.85) 


with  similar  relations  for  the  SV  motion.  In  the  isotropic  limit  ur,  — ►  0  ;  ugp  — ►  0.  Consequently, 


F(oj )  .  e  ,kbR  F(u  e  ,k% 

—  sing — - —  ;  ur„-»  — - — rcos 0 — — 

4*/i  R  Hp  4?r(A  +  2/i)  R 


(4.86) 


II.  Localized  force  in  the  isotropy  plane  (x  direction,  say) 

In  this  case,  the  displacement  held  becomes  more  complicated  due  to  the  presence  of  SH  waves, 
in  addition  to  the  usual  quasi-transverse  and  quasi-longitudinal  waves.  These  are  obtained  directly 
from  (B.35): 


u 


SH 

A 


U 


SH 

<t> 


l£V)  cos<t>  e-xksliR 
drprjsHa5  ksH&2 

tF(w)  ..ft  1 

- -  sin  0  U-  +  - - --j 

IxpVSHQS  l#  *S//A2J 


°  1 


(4.87) 

(4.88) 


where  ksH  is  defined  in  (B.15).  Note  the  line-source  singularity  at  A  =  0.  [This  singularity  arises 
from  the  incongruency  of  the  axial  symmetry  of  the  medium  and  the  geometry  of  the  localized 
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source.  It  is  already  inherent  in  the  operator  V(-2  in  Eq.  (4.23).]  The  P-SV  part  of  the  field  is 
evaluated  from  (4.57)  and  (4.58)  under  the  conditions  ( k,R )  1. 

III.  Explosion  (center  of  compression) 

Since  G  is  symmetric,  displacements  are  given  by  ( Ben-Menahem  and  Singh,  1981) 

u  =  -M0divr  G  (r|r0;u;)  ,  (4.89) 


where  M o  is  the  source’s  moment.  Applying  the  operator  divr  to  (4.43),  we  obtain 

S  =  ~  {^r[w2  +  (a2  -  a3)d]  +  QsV?]'?  +  e2d2[(a3  +  a5  -  a2)d 2  +  (a,  -  a3  -  a5)V2)]'i»|  . 


In  the  isotropic  limit,  the  r  ■  h  ■  s  of  (4.90)  goes  to 


M0l 


U isotropic  =  -^Vr(u/  +  /?2V2)tf  . 

The  two  fields  differ  in  two  respects: 

1.  In  azimuthally-isotropic  media  the  explosion  sources  produce  also  shear  waves. 

2.  The  field  is  no  more  a  gradient  of  a  potential. 

In  the  isotropic  limit  u2,  —  0  u2p  —  -  4nffi2fi)  as  it  should.  Also,  ug ,  - 

ur,  — 1 •  0,  and  ugp  — *  0.  In  the  far-field,  one  obtains  from  (4.90): 


(4.90) 


(4.91) 


0, 


UpA  = 


u,2  - 


U,A  = 


M0ai 

P 

Mpai 

P 

M0a5 

P 


DVr 


Oz  ' 


-91  + 


a  i 


02  ~  «3 
ai 


c92  +  kp 


d^p 

0A  ' 


Q1  ~  Q3a2  ,  02  ,  ,2]  Ws 

— , 4  +  d- +  H  ar- 
01  +  "zzzzoi  +  *2  . 

05  uA 


(4.92) 
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•vhere  k2  =  uj2/ai  ;  k]  =  u>2/as. 


IV.  Torque  (center  of  rotation) 

Let  the  moment  of  the  torque  align  with  the  z-axis.  The  displacements  are  then  given  by 
(Ben-Menahem  and  Singh,  1981) 


u  —  Mq£z  ■  curlT  G  • 


(4.93) 


Considering  (4.38),  it  is  apparent  that  the  term  aaVV'P  does  not  contribute  to  the  field.  Then, 
using  the  tensor  identity  [for  a  fixed  vector  m  and  variable  symmetric  dyadic  G  with  trace  p] 


in  ■  curl  G=  curl(m ■  G)  ~  in  X  div(G  -pi), 


(4.94) 


one  can  show  that  the  rest  of  G  contributes  the  displacement 

i=i^Br,[sradh^)(ksHkUi-]' 

which  renders  SH  waves  only. 

The  radiation  pattern  in  the  absolute  far-field  is 


-iksnMo 


v  4npr}sHas  Lsin2  9  4-  q2  cos2  6}  R 


e-'ksHR 


V.  Strike-slip  displacement  dislocation 


(4.95) 


(4.96) 


The  displacement  field  is  of  the  explicit  form  (Ben-Menahem  and  Singh,  1981) 

«=-A/oM  er.^  +  ey-^  ,  (4.97) 

where  the  slip  is  in  the  a-2-plane  with  a  normal  in  the  y  direction  and  a  moment  Mq(u>).  Applying 
the  operator  (er^  +  ey§^j  to  the  explicit  expression  of  G  in  (4.57)-(4.58),  we  can  derive  the 
various  field  components. 


84 


Note  that  in  all  expressions  for  the  displacements  we  have  not  included  the  contribution  of  the 
singularity  tensors  Go  and  G'o  [Eqs.  (4.51)-(4.52)].  Although  these  terms  are  of  order  j  in  the 
axial  far-field,  their  values  relative  to  the  absolute  far-field  are  less  than  a  few  percent  for  all  the 
crustal  models  listed  in  Thomsen  (1986).  The  use  of  the  integral  representation  (4.80)  enables  one 
to  represent  the  displacement  fields  in  a  form  of  a  Hankel  transform.  The  kernels  of  these  integrals 
for  different  seismic  sources  are  shown  in  Table  1. 


4.5  The  Exact  Total  Displacement  Field 


In  Table  1,  we  have  considered  the  far-field  displacements  for  several  seismic  sources.  In  this 
subsection  we  show  the  exact  total  displacement  field.  In  the  case  of  a  center  of  compression 
(explosion),  we  obtain,  for  the  2  component 


M0{u) 


too 

I  dkh0(z;k)J0(kA), 
Jo 


(4.98) 


where 


M*;*)  =  lzl5^t)-[(a1-a3)A:2-a55?(fc)-tu2]e"|z|5l( 

(4.99) 

For  a  horizontal  force  in  the  z-direction  we  find  for  the  A  component  of  the  P-SV  displacement 


(4.100) 


whore 


hi{z;k)  = 


VAk*  +  B  F  +  E 2 


[o5fc2  -  a25|(fc)  -  a;2]- 


-MS2(fc)  „ 


-|*|S,(fc) 


The  expression  for  u&  in  (4.100)  must  be  supplemented  by  the  corresponding  SH  counterpart  (4.87). 

The  fields  of  other  sources  can  be  represented  in  a  similar  way.  When  this  source  representation 
is  coupled  to  a  multilayered  structure  by  means  of  Haskell-type  algorithms,  one  can  compute  the 
fields  of  point  sources  in  inhomogeneous  anisotropic  media. 


4.6  Extension  of  the  Theory  to  Sources  in  Multilayered  Media 


In  multilayered  isotropic  media,  the  spectral  response  to  point  sources  is  constructed  from  the  known 
solution  of  a  homogeneous  layer  and  the  source  integral  representation  in  a  boundless  domain  [e.g., 
Haskell,  1964],  We  wish  to  apply  this  matrix  algorithm  to  azimuthally  isotropic  media  in  which  the 
axis  of  symmetry  coincides  with  the  vertical  axis.  Let  us  therefore  write  the  fundamental  potential 
of  the  /- th  layer  in  its  Ilankel  transform  representation 


1  f°° 

*(°  =  - 777-m  /  h{l)(z;k)J0(kA)kdk  ,  (4.102) 

8x02 'ag  ^0 

where  h^l^(z,k)  is  given  explicitly  in  (4.80).  We  then  substitute  (4.102),  together  with  the  known 
representation 


roo 

6(f—  f0)  =  6{z  -  z0 )  /  J<j(kA)kdk  , 
Jo 

in  (4.36).  The  resulting  equation  for  h^(z-,k)  is  then 


(4.103) 


-o) 


(1-104) 
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where  the  {d^\  d \j\  d^}  are  defined,  for  each  layer,  in  (4.79).  To  obtain  the  surface  displacements 
for  either  surface  waves  or  body  waves  on  the  surface  of  a  multilayered  half-space,  one  follows  the 
procedure  outlined  in  Haskell  (1964)  and  Ben-Menahem  and  Harkrider  (1964).  To  this  end,  we  need 
the  explicit  general  solution  (4.104)  at  z  ^  zq  in  the  form  of  a  sum  of  four  wave-types  corresponding 
to  upgoing  and  downgoing  quasi-longitudinal  and  quasi-transverse  plane-waves.  To  these  we  add 
the  source  representation  in  the  source  layer  given  via  (4.102)  for  the  particular  source  in  question. 
One  can  then  calculate  the  fields  of  point  sources  in  a  multilayered  anisotropic  media. 
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5.  STATIONARY  PHASE  APPROXIMATION  OF  THE  GREEN’S  TENSOR 


We  shall  next  derive  the  asymptotic  Green’s  tensor  for  lot  — *■  oo.  When  ur  is  large,  the  contributions 
from  the  integration  over  complex  angles  (inhomogeneous  waves)  are  negligible  compared  to  the 
contributions  from  the  homogeneous  waves.  Therefore,  the  important  contributions  to  the  integral 
in  (3.20)  will  arise  from  points  on  the  real  surface  (Buchwald,  1959;  Lighthill,  1960). 

To  this  end  we  use  the  method  of  stationary  phase,  which  says  that  (3.20)  is  asymptotic  as 
ur  — ►  oo  to  a  sum  of  contribution:  from  points  on  the  slowness  surface  l?m  where  r  ■  s  is  stationary, 
namely,  from  points  on  the  slowness  surface  where  the  normal  is  parallel  to  r. 

Denoting  these  points  as 

sO),  i<2),  i<3>,  ...,  s<">  ,  (5.1) 


we  obtain  the  asymptotic  solution 


where 


C(^l))  =  -i  exp  j-^[s£ra  (*!(#))  +spn  («2(«<')))]|  •  ( 

Here  «i(sW)  and  k2 (^)  are  the  principal  curvatures  of  tne  slowness  surface  at  s  =  i^),  and 


=  «l(5<,))K2(5<,))  , 


(5.4) 


is  the  Gaussian  curvature  of  the  slowness  surface,  assumed  non-zero,  evaluated  at  s  =  The 
principal  curvatures  kx  and  rc2  are  positive  if  the  surface  is  concave  to  the  direction  r.  The  case  of 
parabolic  points  of  the  slowness  surface  at  which  the  Gaussian  curvature  is  zero  can  also  be  handled 
using  the  stationary  phase  approximation,  but  will  not  be  considered  in  this  paper. 


We  can  see  from  (5.2)  that  the  asymptotic  expression  of  the  Green’s  tensor  depends  on  the  group 
velocity  and  the  curvature  of  the  slowness  surface.  Also,  we  note  that  the  amplitudes  are  inversely 
proportional  to  the  distance  from  the  source  for  all  waves  corresponding  to  ordinary  points  of  the 
slowness  surface.  The  error  in  (5.2),  is  £>(l/r2)  .  Besides,  we  have  assumed  that  the  integrand 
in  (3.20)  varies  slowly  near  the  stationary  points.  The  availability  of  the  Green’s  tensor,  in  its 
exact  (equation  (3.20))  and  asymptotic  form  (equation  (5.2))  enables  one  to  obtain  analytical  and 
numerical  solutions  for  a  wide  range  of  media  and  sources. 

Location  of  Stationary  Points 

We  now  set  to  evaluate  the  stationary  points.  Let 

P(s)  =  jA  =  r-5 ,  (5.5) 

m 

in  the  argument  of  the  exponential  in  (3.20),  anu  let 

rm(s)  =  0  ,  (5.6) 

be  the  slowness  surface  equation.  Since  the  stationary  points  of  the  integral  in  (3.20)  must  belong 
to  the  slowness  surface,  it  follows  that  the  stationary  points  satisfy 

V;[p(s)  -  t'rm(5)]  =  0  ,  (5.7) 

Mi)  -  0  ,  (5.8) 

where  v  is  a  Lagrangian  multiplier.  Since 

V?p(s*)=f,  (5.9) 
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and 
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with  s  =  (s2  +  s \)x!2. 

Due  to  the  azimuthal  symmetry,  we  consider  f  =  A£a  +  zez  and  the  stationary  points  are 
obtained  by  solving  the  system 


Thus,  we  find 


S3 


A  -  v 

2? 
»  3 

II 

o 

z  -  V 

=  0 
ds3 

(5.17) 

r 

O 

II 

£ 

{~T 

\C6eJ 

sin# 

(5.18) 

(sin2  0  +  ^  cos2  o')  ^ 

cos  0 

C  44 

(5.19) 

/  \  1  /2  ’ 

(sin2  0  +  cos2  0) 

where  6  is  the  collatitudinal  angle  corresponding  to  the  observation  point. 


P-SV  case: 

The  slowness  surface  is 


rTO(.s)  =  m  +  eVU  -  2  , 


(5.20) 


where  e  =  1  in  the  P  case  and  e  =  -1  in  the  SV  case,  and 


M  = 


m=  (fulfil)  ,2+  s23i 

^Cn  -  C44j  a2+  ^33  ~C44j  s2j2_4^£llJl£44^  ^33  -  C44 


(5.21) 


)-( 


c13  +  c44\  2 

p  ) 


„2„2 

s  s3 


(5.22) 
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The  stationary  points  are  obtained  by  solving  the  system  (5.17)  with  rm(s)  given  in  (5.20). 


Eliminating  z/,  we  have  the  system 

.  r  dm 


+ 


dM  i 


[dsz  2 n/M  ds3 


f  dm  € 


dM 


Z  [  ds  ^  2 %/M  9s 

rm(iO  =  o . 


=  o, 


(5.23) 


The  exact  solution  of  (5.23)  for  the  stationary  points  can  only  be  achieved  numerically.  This  task 
is  now  under  study.  In  order  to  obtain  analytical  solutions  for  the  stationary  points  we  consider  an 
approximate  form  of  the  slowness  surface  (or,  equivalently,  of  the  phase  velocity),  the  justification 
of  which  is  stated  in  sections  5.2  and  5.3.  We  assume 

i  \ 

(5.24) 


=Cf(s2  +  '  1  • 


in  the  P  case  and  in  the  SV  case 


Tsv(s)  =  ~  (*2  +  ~  1  ’  (5.25) 

Here,  tjp  and  r),  are  given  later  in  (5.42)  and  5.60)  respectively.  Using  these  forms  of  the  slowness 
surfaces,  we  find  for  the  stationary  points: 

Quasi-longitudinal  wave: 


,  _  .  (±\  sing 

\cu  )  (sin2  9  +  vj  cos2  0)1/2  ’ 

f_p_\  tjIcosO 

53  \cn  /  (sin20  -f-  tfcos2#)1/2 

Quasi-transverse  wave: 

(p  \  sin# 

C44  )  (sin2  0  +  rfi  cos2  By!2  ' 
.  ( _P_\  r,]  cos  9 

53  VC44  /  (sin2  6  +  tj]  cos2  0)1/2 


(5.26) 

(5.27) 


(5.28) 

(5.29) 
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The  particular  case  of  ellipsoidal  anisotropy  (Gassman,  1964)  can  be  easily  obtained  from  the 


previous  results  by  setting 


2  *"tl  j  ■, 

Vv=  —  and  T),  =  1, 


(5.30) 


in  (5.24)-(5.29). 


5.1  SH  Motion 


We  shall  estimate  the  asymptotic  SH  Green’s  tensor,  using  the  results  obtained  in  Section  5.  Since 
the  slowness  surface  is  an  ellipsoid,  there  is  only  one  point  on  it  where  the  normal  is  parallel  to 
a  given  direction.  Therefore,  if  ( x,y,z )  are  the  coordinates  of  a  given  observation  point  in  space, 
there  will  be  only  one  kind  of  wave  passing  through  it,  corresponding  to  a  point  (s^, s^)  on 
the  slowness  surface.  The  coordinates  of  the  stationary  point  that  give  rise  to  outgoing  waves  are 
obtained  from  (5.18)  and  (5.19)  by  choosing  the  positive  sign. 

The  evaluation  of  the  Gaussian  curvature  is  simplified  when  the  slowness  surface  is  axially 
symmetric.  It  is  sufficient  to  evaluate  1C 3  in  the  case  s 2  =  0,  and  owing  to  the  symmetry,  the  value 


thus  obtained  holds  for  all  values  of  $2-  We  obtain 


IC3(?)  = 


Evaluating  Eq.  (5.31)  at  the  stationary  point  v»e  *>ave 


£  («?  +  *2  +  ^5) 


(5.31) 


£3(5^)  =  (sin20  + —cos2#') 
P  \  c 44  / 


(5.32) 


The  magnitude  of  the  group  velocity  at  the  stationary  point  can  be  obtained  from  (4.5).  We 
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have 


w3(^l))  =  ^  ^sin20  +  — cos2#^ 


-1/2 


(5.33) 


«(3) 

The  projection  tensor  A  can  be  written  conveniently  in  cylindrical  coordinates  as 


~(3)  ,  -s 

A  («)  =  . 


(5.34) 


In  this  case  both  principal  curvatures  are  negative  so  that  C(s^)  =  1.  Therefore,  the  asymptotic 
expression  of  the  Green’s  tensor  for  SH  waves  becomes 

2  exp  to >(/>/c66)1/2  (sin20  +  ^cos2^12  rj 

MC44C66)1/2  (sin2<,+  ^co52^1/2r 

or,  since  R  =  r  ( sin26  +  l^-cos2#)  ^  ,  (see  (4.26)),  we  can  write 


<-,SH 

G  (r,w)  = 


-e4>e<t> 


(5.35) 


~sh  1 

G  (r,u>)  = 


exp  [~tu;(p/c66),/2^]  . 


&4>^4>  —  ff(^t  A,  u>)e^,e^  ,  (5.36) 


4ir(c44c66)»/2  J? 

where  we  have  used  the  definition  of  5(2,  A)  given  in  (4.25).  Using  (4.93)  we  have  in  homogeneous 
anisotropic  media 


-Mo 


rffrocfr 


B-«w(p/cee)1/J>-' 


x  e* 


’r?ue  ~  47r(c44c66)i/2 
Therefore,  in  the  far  field  we  have  for  the  displacement  field  of  SH  waves 


(5.37) 


,  -iw(p/c66)1/2M0  | 

f  sind  ^ 

exp 

1 

-iu(p/c66)1/2  1 

^ sin2d  +  j^-cos2^j 

1/2 

1  r 

tor,ue  4*(c44c66)V2  1 

^sin2#  +  sjs. cos20  J 

(5.38) 


in  agreement  with  Ben-Menahem  (1989). 
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5.2  Quasi-Longitudinal  Motion 

The  asymptotic  expression  corresponding  to  the  P  part  of  the  total  Green’s  tensor  (equation  (3.21)) 


Gr  (f»  =  J-£  ..MnL  it'” 

4x pr  “  u)p(st'})|^p(it'))|i/2 


(5.39) 


We  have  stated  earlier  that  the  analytic  determination  of  the  stationary  points  in  the  case  of 
azimuthally-isotropic  media  is  hampered  by  the  special  algebraic  form  of  C2  (equation  (4.66)). 

We  now  consider  an  approximation  of  the  phase  velocity  that  enables  us  to  find  closed  form 
expressions  for  the  stationary  points.  Assume  that  we  can  write  the  P  phase  velocity  as 


Cl  =  a\ (sin2  9  +  -=•  cos2  6)  , 


(5.40) 


where  q2  is  still  unknown.  It  then  follows  algebraically  from  (4.66)  that 


\/~M.  —  (ai  -  05)  sin2  9  +  ~  a2  ~  a5 J  cos2  9  . 


(5.41) 


When  we  compare  (5.41)  with  M  in  (4.66)  without  the  last  term  4a 
the  solution  for  the  unknown  r?2  is 


r^(T3;:^r-a5)cos4i 


2  _  Qi(Qi  -  as) 
p  05(01  -  05)  +  a23  ' 


(5.42) 


We  then  conclude  that  our  approximation  is  valid  whenever  the  neglect  of  the  above  term  in 
the  expression  of  the  phase  velocity  is  justified.  For  isotropic  and  ellipsoidally  anisotropic  media 
(equation  (4.67)),  this  term  is  identically  zero.  In  Figure  3a,  we  show  the  relative  error  in  the  quasi¬ 
longitudinal  phase  velocity  when  using  this  approximation  for  some  earth  materials  considered  by 
Thomsen  (1986)  (see  Table  2).  We  can  see  that  this  approximation  is  uniformly  valid  for  a  wide 


range  of  angles  and  the  error  increases  as  we  approach  the  axis  0  =  0. 
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i,From  (5.40)  we  obtain  the  equation  of  the  slowness  surface  for  quasi-longitudinal  waves 


ai(si  +  32  +  ~2S3)  -  1  • 
vp 


(5.43) 


This  surface  is  an  ellipsoid,  and  in  the  case  of  isotropy  (7?p  =  1),  it  reduces  to  a  sphere  of  radius 


_L  _  p 

ai  \j  A+2/i  ’ 


The  coordinates  of  the  stationary  point  that  correspond  to  outgoing  waves,  are  obtained  from 
(5.26)  and  (5.27)  by  choosing  the  positive  sign.  In  Figure  4,  we  show  the  geometry  of  the  slowness 
surface  in  a  generic  case  for  which  r)v  >  1. 

The  magnitude  of  the  group  velocity  is  easily  found  by  differentiation  of  the  expression  of  the 


phase  velocity 


wp(s)  =  ai  h  +  — 

\  Vp 


2  .  1  „2 


(5.44) 


Evaluating  (5.44)  at  the  stationary  point,  we  obtain 


(sin2  9  -f  7?2  ros2  9)1!2 


(5.45) 


Next,  we  calculate  the  Gaussian  curvature.  The  axial  symmetry  simplifies  the  evaluation  and  we 


At  the  stationary  point 


1  s2  +  si  In2 


/Cpf^1')  =  ^-(sin2  6  +  r)l  cos2  0 )2  . 


(5.46) 


We  consider  the  eigenvector  corresponding  to  quasi-longitudinal  waves.  For  azimuthal  isotropy. 


we  can  write,  after  some  algebra,  the  displacement  polarization  eigenvector  at  the  stationary  point 
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=  sin  ape&  +  cos  apez  , 

(5.48) 

.  MO) 

s,nvp  =  s >*»m, 

(5.49) 

2  MO) 

cos  op  -  Ti*cosO  , 

(5.50) 

di(6)  =  (ax  -  a5)  sin2  9  +  [ax  -  772(a2  -  a3)>72  cos2  9  , 

(5.51) 

d2(9)  =  a3  sin2  9  +  (ax  -  a^rj^rjl  cos2  9  , 

(5.52) 

D(9)  =  [d2(<?)  sin2  9  +  d\{9)T]*  cos2  0]1/2  . 

(5.53) 

In  the  isotropic  limit  crp  =  9,  and  we  obtain  the  usual  form  of  the  displacement  polarization  vector 
for  P  waves. 

Using  (5.43)  we  have  for  the  projection  operator 


‘-(p) 

.4 


=  =  sin2<rpeAeA  +  sin<7pcosc7p(e^ez  +  eze&)  +  cos2  crpezez 


(5.54) 


Then,  substituting  (5.45),  (5.47)  and  (5.54)  into  (5.39),  and  using  the  fact  that  C(i^)  =  1,  we 
obtain  the  asymptotic  expression  of  the  Green’s  tensor  for  quasi-longitudinal  waves 


Gp  (r,v)  = 


rjp  ~(p)  e 
47rpai  A  Rp 


(5.55) 


where 

Rp  =  r(sin2  9  +  r)2  cos2  9)1/2  . 


(5.56) 
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5.3  Quasi-Transverse  Motion 


In  this  case  the  asymptotic  expression  of  the  Green’s  tensor  is 


1  N 

Gsv  (r,u)  =  - - 

A  ir  tvr 


C(W) 


Airpr  o;5v(iI9)|£sv,(stO)Jl/2 


A  (s<l))H(r- ^)e 


—  i 


(5.57) 


and  the  corresponding  phase  velocity 


Csv  -  7t(m  ~  v  M)  i 


(5.58) 


with  m  and  M  given  in  (4.66). 

At  this  point  we  follow  the  same  arguments  considered  for  approximating  the  phase  velocity  in 
the  quasi-longitudinal  case.  Assuming  that  we  can  write  the  SV  phase  velocity  as 


Csv  =  a5(sin2  9  +  —  cos 2  9)  , 

n* 


(5.59) 


with  q2  still  unknown,  it  follows  algebraically  from  (5.58)  that 


2  _  <15(01  -  Q5) 

3  02(01  -  05)  -  a| 


(5.60) 


as  long  as  the  last  term  in  M  in  (4.66)  can  be  neglected.  In  Figure  3b  we  show  the  relative  error  in 
phase  velocity  when  using  the  approximate  expression  (5.59)  for  some  earth  materials  considered 
by  Thomsen  (1986)  (see  Table  2).  Again,  we  observe  that  the  approximation  is  uniformly  valid  for 
a  wide  range  of  angles  except  along  the  axis  9  =  0. 

^From  (5.60)  we  know  that  this  approximation  for  the  quasi-transverse  waves  fails  whenever 


02(01  -  o5)  <  a2  , 
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yielding  unbounded  or  imaginary  values  of  t)3  at  those  points.  Although  (5.61)  is  not  excluded  by 
the  physical  realizability  conditions,2  media  of  this  type  are  rare  in  earth  sedimentary  rocks  (only 
2  out  of  58  cases;  Thomsen,  1986). 

The  equation  of  the  slowness  surface  is  obtained  from  (5.59) 


a5  I  +  $2  +  1  —  1  > 


(5.62) 


with  the  corresponding  stationary  point  for  outgoing  waves  given  in  (5.28)  and  (5.29),  choosing 
again  the  positive  sign.  The  magnitude  of  the  group  velocity  at  the  stationary  point  takes  the  form 


wsv(s{1))  = 


(sin2  6  +  rfi  cos2  0)1/2 


(5.63) 


and  the  Gaussian  curvature  is 


ICsv(^)  =  -f(sin2  9  +  rjl  cos2  9)2 


(5.64) 


Next,  we  consider  the  displacement  polarization  eigenvector  for  the  quasi-transverse  waves  at  the 
stationary  point.  After  some  algebra,  one  obtains 


rrSV 

1  =  cos o,e&  —  sin<7,ez  , 


(5.65) 


where 


cos  crs  =  r}a  cos  0 


E(9 )  ’ 


E{9)  ' 


(5.66) 


(5.67) 


2CuC33  >  Ci  3,  C  ll  >  C33,  (cil  +  Cl  2  )C33  >  2c?3 
c  11  >  0.  c 33  >  0,  c, 4  >  0.  cgc  >  0,  c  1 1  >  (r  1 2 1 
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(fii  -  as)  sin2  9  +  | j^-j  a1-q2(a2-a3)  q]  cos2  9 


=  a3sin2  9  +  \^~ J  ai  “  asV*  q]cos29, 

=  [62(0)  sin2  9  +  b2(9)q 3  cos2  9 j1^2  . 


(5.68) 


(5.69) 


(5.70) 


In  the  isotropic  limit  crs  =  9  and  the  conventional  expression  of  the  displacement  polarization  vector 
for  SV  waves  follows. 

i,From  (5.65)  the  corresponding  projector  operator  is 

A  =  f(sv)f(S\  )  _  cos2  aae&e&  -  sin<r3  cos  a,{e^ez  +  eze&)  +  sin2cr3e2eJ  .  (5.71) 

Therefore,  using  (5.63),  (5.64),  (5.71)  and  (5.57),  together  with  the  fact  that  C(s^^)  =  1,  the 
asymptotic  Green’s  tensor  for  the  quasi -transverse  SV  waves  is 


«  n2  ~(SV)  g-xwR./aJ 

=  wsr 


(5.72) 


where 


Rs  =  r(sin2  9  +  rj2  cos2  9 )1^2 


(5.73) 


5.4  Divergence  Coefficients  and  Wave  Surfaces 


We  have  noticed  in  (5.36),  (5.55),  and  (5. 72)  a  typical  divergence  coefficient  for  the  Green's  tensor 
of  the  form 


r(sin2  9  +  T)2  cos2  9 )1/'2 


(5.74) 


where  q  =  qp  qs,  or  qsii-  In  order  to  give  a  geometrical  interpretation  to  the  angular  dependence 


of  (5.74),  we  consider  the  corresponding  wave  surface.  Its  equation  can  be  obtained  by  following 
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waves  for  which  the  phase  u>(r  ■  s  —  t)  is  zero,  then,  at  t  —  1, 


r  ■  s  =  1  . 


(5.75) 


The  waves  at  point  fin  space  are  generated  by  a  point  s  of  the  slowness  surface,  satisfying  Tm(f)  = 
0,  at  which  the  normal  is  parallel  to  the  vector  f,  namely, 


f  =  qv rm  , 


(5.76) 


q  being  a  scalar.  Taking  the  scalar  product  of  (5.76)  with  s  and  using  (5.75),  we  obtain 


VTm  •  s 


(5.77) 


Therefore,  the  coordinates  of  a  point  of  the  wave  surface,  at  t  =  1,  will  be  given  by 


vrm  w 

r  =  — - =  — — -  =  w  , 

vrms  W-S 


(5.78) 


where  w  is  the  group  velocity  vector. 

Since,  under  the  approximations  considered  in  this  paper,  the  P,  SV  and  SH  slowness  surfaces 
are  ellipsoids,  it  is  sufficient  to  exhibit  only  one  example.  For  the  P  case  it  follows  from  (5.78)  that 
the  equation  of  the  wave  surface  is 


— (x  +  jT  +  )  =  1  • 


(5.79) 


In  Figure  5  we  show  the  geometry  of  the  wave  surface  in  a  generic  case  for  which  qp  >  1.  We  can 


easily  see  that 


a  /2  ^isotropic  R 

T=T  =  V-,-  =  sin2  9  +  1)1  cos2  0)’/2  =  -P 

uJ(p)(s)  lp  '  r 
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Therefore,  we  find  that  the  ratio 


iposition  on  isotropic  w_ave_surfaccL  Q  ivalcnUy,  the  ratio 
Iposition  on  anisotropic  wave  surfacel  n 


[magnitude  of  isotropic  group  velocity! 
[magnitude  of  anisotropic  group  velocuyi’ 


along  the  direction  of  the  observation  point  gives  the 


angular  dependence  of  the  divergence  coefficient  (5.74). 


5.5  Anisotropy  Data 

In  order  to  apply  our  theory  we  consider  an  exhaustive  listing  of  measurements  of  anisotropy  for 
58  earth  materials,  including  mostly  sedimentary  rocks,  reported  in  Thomsen  (1986).  These  data 
were  interpreted  by  the  original  investigators  in  terms  of  the  five  elastic  constants  of  azimuthal 
isotropy.  About  70%  these  models  have  a  maximum  error  below  30%  at  9  —  0  and  180°when  we 
use  the  approximate  expression  for  the  phase  velocities  given  in  (5.40)  and  (5.59).  However,  in  the 
range  30° <  6  <  150°,  the  error  is  approximately  below  10%  (see  Figure  3).  In  Table  2  we  show  a 
partial  list  of  these  models,  including  the  elastic  constants,  77*,  77^,  and  the  anisotropy  parameter  c 
(percentage) 

e(%)  =  -n.-C33  x  100  ,  (5.81) 

■^33 

as  a,  somewhat  arbitrary,  measure  of  anisotropy. 

As  another  indicator  of  the  fit  of  the  approximation  of  the  phase  velocities,  we  show  in  Figure  6 
the  exact  and  approximate  slowness  surfaces  for  four  selected  models  (see  Table  2).  The  good  fit 
of  the  approximation  for  a  wide  range  of  angles  centered  at  0  =  90°is  again  apparent.  Besides,  we 
can  see  that  this  approximation  induces  ellipsoidal  slowness  surfaces,  and  therefore  excludes  the 
possibility  of  parabolic  points. 

In  the  next  subsection  we  evaluate  the  displacement  fields  and  radiation  patterns  of  some 
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common  seismic  sources. 

5.6  Radiation  Patterns 

I.  Localized  force  in  </u  direction  of  the  symmetry  axis. 

The  representation  of  this  field  in  terms  of  the  fundamental  potential  $  was  discussed  already  in 
Section  4.3.  It  can  also  be  given  directly  in  terms  of  the  Green’s  tensor 

u  =  F(u)[Gp  +  Gsv]  •  e2  .  (5.82) 

Using  the  asymptotic  expressions  in  (5.55)  and  (5.72),  we  obtain  the  results  shown  in  Tables  3  and 

4. 

In  Figures  7a  and  7b  we  show  the  vertical  radiation  patterns  of  the  amplitude  spectrum,  namely, 
the  dependence  of  the  spectral  displacement  field  of  uzp  and  u23 ,  respectively,  on  6  for  fixed  azimuth 
and  r.  The  isotropic  limit  is  marked  by  a  broken  line.  The  deviation  of  the  results  from  the 
isotropic  case  demonstrates  that  these  patterns  can  be  used  as  diagnostic  aids  for  the  identification 
of  anisotropic  regions  in  the  earth’s  crust  and  mantle. 

II.  Localized  force  in  the  isotropy  plane  (x- direction). 

Here,  the  displacement  field  is  determined  by 

F(u)[GP  +  Gsv]-ec  ,  (5.83) 

and  using  (5.55)  and  (5.72)  we  obtain  the  results  shown  in  Tables  3  and  4. 

In  Figures  8a  and  8b  we  show  the  vertical  radiation  patterns  corresponding  to  the  amplitude 
spectrum  of  ugp  and  ur3,  respectively.  These  components  are  null  in  the  isotropic  limit,  and  may 
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thus  be  used  to  diagnose  azimuthal  isotropy  provided  a  powerful  enough  source  is  used  and  the 
collatitudinal  coverage  in  the  recordings  is  adequate. 

III.  Explosion  (center  of  compression) 

Since  G  is  symmetric,  the  displacement  field  is  given  by  (4.92).  Using  the  asymptotic  expressions 
for  the  Green’s  tensor  in  (5.55)  and  (5.72),  we  find  the  results  indicated  in  Tables  3  and  4. 

Contrary  to  the  isotropic  case,  explosion  sources  in  azimuthally  isotropic  media  also  produce 
shear  waves.  In  Figures  9a  and  9b  we  show  the  vertical  radiation  patterns  corresponding  to  the  uzp 
and  u23  (not  present  in  the  isotropic  case)  components,  respectively.  Using  these  two  components, 
we  can  evaluate  the  SV/P  ratio  of  the  displacement  field.  The  results  are  shown  in  Figure  10.  We 
observe  that  SV/P  ratios  as  large  as  4.5  are  predicted  (Model  5,  Robertson  and  Corrigan,  1983). 
This  indicates  that  explosion  sources  can  be  used  diagnostically  to  establish  anisotropy  provided 
that  sufficiently  strong  quasi-longitudinal  waves  are  generated. 

In  Figures  11a,  lib,  12a,  and  12b,  we  present  the  vertical  radiation  patterns  corresponding  to 
the  ua3,  u$p,  lie,  and  uTt  components,  respectively.  All  these  components  are  null  in  the  isotropic 
limit.  For  anisotropic  media  we  have  a  variety  of  new  radiation  patterns  including  an  “eight  lobe” 
pattern  for  Model  3  (Lin,  1985).  This  behavior  can  therefore  be  used  effectively  in  the  identification 
of  azimuthal  isotropy. 

IV.  Strike-slip  displacement  dislocation 

The  displacement  field  in  this  case  has  already  been  given  in  Section  4.3.  Applying  the  operator 
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(e»  ■  +  ey  ■  to  the  asymptotic  expressions  (5.55)  and  (5.72),  we  derive  the  various  field 

components.  The  corresponding  results  are  given  in  Tables  3  and  4. 

In  Figure  13  we  compare  the  radiation  patterns  of  ugp  for  two  sources  in  tuff  (Model  7,  Table  2). 
The  analytic  form  of  this  pattern  is  listed  in  Table  4.  The  sources  are:  a  point  strike-slip  dislocation 
in  the  xz  plane  and  a  point-explosion.  The  patterns  are  almost  identical  in  shape. 

In  Figure  14  we  compare  the  radiation  patterns  of  the  same  sources  in  Mesa verde  clayshale 
(Model  1,  Table  2).  Here  we  display  the  u2a  displacement  component  (Table  3).  Again,  the 
patterns  are  indistinguishable. 

These  results  show  that  even  in  weak  anisotropic  media,  explosions  can  masquerade  for  strike- 
slip  earthquake  fauits.  The  ramifications  of  these  results  are  discussed  elsewhere  (Ben-Menahem 
and  Toksoz,  1989). 

Note  that  our  radiation  patterns  were  calculated  for  the  case  where  both  source  and  receivers 
are  in  the  same  anisotropic  medium.  This  may  seem  a  rather  specific  case.  Nevertheless,  the 
results  are  useful,  among  other  applications,  in  geophysical  exploration  for  oil.  For  example,  in  the 
situation  where  the  source  is  located  in  an  open  hole  section  of  a  borehole. 
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6.  DISCUSSION  AND  CONCLUSIONS 


We  have  derived  analytic  approximations  for  the  far  field  of  point  sources  in  azimuthally-isotropic 
media.  Our  results  exhibit  the  strong  effects  that  anisotropy  has  on  the  source’s  directional  energy 
radiation  and  the  partition  of  energy  between  the  longitudinal  and  shear  components  of  the  motion. 
We  have  also  shown  that  Haskell  matrix  method  is  applicable,  without  modification,  to  multilayered 
azimuthally  isotropic  media  in  which  the  axis  of  symmetry  and  inhomogeneity  coincide. 

In  our  study  we  have  employed  the  multidimensional  stationary  phase  approximation  which 
limits  the  results  to  frequency  and  spatial  domains  such  that  u >r  is  large.  Also,  in  the  process  of 
determining  the  stationary  points  we  found  it  convenient  to  approximate  the  phase  velocities  itself. 
The  price  payed  for  this  artifice  was  the  lesser  accuracy  of  the  spectral  amplitudes  on  and  near  the 
axis  of  symmetry.  This  restriction  can  be  removed  by  solving  numerically  for  the  stationary  points. 
The  alternative  for  the  stationary  phase  approximation  is  the  numerical  calculation  of  the  Hankel 
transform  in  the  exact  analytical  solution.  Note  that  by  our  simultaneous  use  of  the  stationary 
phase  approximation  and  the  approximated  phase  velocities  we  avoid  the  parabolic  points  on  the 
slowness  surface  since  these  surfaces  are  approximated  by  ellipsoids.  In  general,  however,  the 
stationary  phase  approximation  fails  at  parabolic  points  of  the  slowness  surface  and  a  higher  order 
approximation  must  be  made  at  these  points. 

The  availability  of  the  Green’s  tensor  in  closed  form  enables  one  to  extend  the  results  to 
anisotropic  half-space  configurations  where  the  axis  of  symmetry  is  not  necessarily  normal  to  the 
free  surface.  This  will  be  used  to  model  azimuthal  anisotropy 
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APPENDIX  A.  DIFFERENTIAL  GEOMETRY  OF  THE  SLOWNESS  SURFACF 

Some  differentio-geometrical  properties  of  the  slowness  surface  are  needed  for  the  evaluation  of  the 
Green’s  tensor.  Let  the  position  vector  of  a  point  on  the  slowness  surface  be  parameterized  as 

si  =  R(0,<j>)sin0  cos<f>  , 

s2  =  R(0, 4>)  sin  0  sin <j> ,  (A.l) 

S3  =  #(#,<£)  cos  0 

with  0  <  0  <  x,  0  <  </>  <  2  x. 

^From  (3.12)  we  have 

\S\  =  pj-  (A.2) 

'-'m 

where  Cm  is  the  phase  velocity  of  type-m  waves.  It  then  follows  from  A.l  and  A.2  that 


W,4>)  = 


Cm(0,<t>)  • 


We  shall  omit  the  subscript  m  in  R(0,<f>)  for  convenience  but  it  is  understood  that  the  function 
R(0,<p)  will  be  different  depending  on  the  sheet  of  the  slowness  surface  under  discussion. 

We  now  evaluate  the  two  vectors  and  which  lie  in  the  tangent  plane.  We  obtain 

80  ~  80er  +  Ree  ’  (A'4) 


8 &  8R 

The  slowness  surface  normal  is  then  obtained  from  (Struik,  1961) 


8s  8s  (  .  8R  .  ..  8R\ 
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;From  A. 6,  the  surface  element  on  the  slowness  surface  can  be  easily  evaluated  (see  (3.17)) 


=*>  da  =  R2  <  sin2  6 


da  =  \N\d6d4> 

I 


1  + 


,  1  dR\2V/2  .... 

+  fI^U  Mi*. 


(A. 7) 


,R  89  ) 

Another  important  quantity  in  evaluating  the  Green’s  tensor  is  the  group  velocity.  We  obtain 
it  by  differentiation  of  the  phase  velocity  of  type-m  waves 

^  .  8Cm 


WK 


86 


Its  magnitude  takes  the  form 

wm{9,4>) 

or.  in  terms  of  R(9,4>), 


C2m  + 


dCm\ 2 
86  J 


“'m(M)  =  Jp 


,/From  (A. 6)  and  (A. 8),  we  also  find 


1 

8Cm . 

^  sin# 

8<t>  €* 

1 

( 8Cm > 

sin2  9 

\  8<j>  ) 

1 

+  .  o 

(A. 8) 


1/2 


(A. 9) 


1/2 


(A. 10) 


JV  = 

^  m 


(A.  11) 


as  expected,  since  the  group  velocity  vector  is  normal  to  the  slowness  surface. 

Finally,  we  give  the  necessary  expressions  for  the  evaluation  of  the  Gaussian  curvature.  We 
start  with  the  first  and  second  fundamental  forms  of  the  slowness  surface  (Struik,  1961) 


I  =  E(d9)2  +  2Fd6d<t>  +  G{d<t>)2  =  ds-ds,  (A. 12) 

II  =  e(d0)2  +  2 fd9dd>+ g(do)2  =  -ds-  X  .  (A.  13) 
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where 


Os  ds 


86  ‘ 

86 

ds 

ds 

86  ‘ 

8<p 

ds 

ds 

8<t> 

84> 

82s 
86 2 

■  N 

82s 

n 

868<j> 

82s 

■1 V 

(A. 14) 
(A. 15) 
(A. 16) 
(A. 17) 
(A.  18) 
(A. 19) 


The  Gaussian  curvature  is  then  given  by 


eg  -  f2 


(A. 20) 


EG-F 2  '  ' 

The  second  derivatives  needed  in  ( A.17)~( A.  19)  can  be  evaluated  using  our  general  parameter¬ 


ization 


82s 

~8(P 

82s 

868o 

82s 

d<t>2 


868<b 

(d2R 

\  8<t> 2 


8R. 

86  €e  • 

(A. 21) 

fdR  n  A 

[-QQ  sind  +  R  cos  6J  , 

(A. 22) 

„  OR  .  n 

er  -  R  sin  6  cos  deg  -f  2 - sin  6c0  . 

do 

(A.  23) 

Thus,  we  have 


(8R\2  n2 

\w)  +R  • 

or  cm 

oe 


(A. 24) 
(A. 25) 
[A. 26) 
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e  = 


/  = 


S  = 


R  sin  9 

|a7I 

R  sin  9 

~w 

R  sin  6 

~w 


■©-»  )-=(»)■ 


d<p  )  09  0<p 


r(UL 

\09d4> 

nl  d2R  _  .  2a  OR  .  -  \  ,  /dtf'2 

<,+  _slnScosej-2(_ 


(A. 27) 
(A. 28) 
(A. 29) 


Substituting  ( A.24)-(A.29)  into  (A. 20),  the  Gaussian  curvature  of  the  slowness  surface  can  be 
evaluated  for  the  general  case. 

The  results  in  this  appendix  can  be  used  to  evaluate  the  Green’s  tensor  in  the  case  of  general 
anisotropy  provided  a  suitable  expression  of  the  phase  velocity  is  obtained. 
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APPENDIX  B.  AUXILIARY  POTENTIALS  FOR  AZIMUTHALLY-ISOTROPIC  MEDIA 


The  Fourier-transformed  Cauchy  equation  of  motion  for  an  elastic  solid  under  a  body-force 


distribution  F(r,u>),  is 


div  T  +pu2u  =  —pF  ; 


where  p  is  the  density,  ui  is  the  angular  frequency,  u(f,u)  is  the  displacement  at  field  point  r  and 
T  is  the  dynamic  stress-tensor.  Choosing  a  cartesian  coordinate  system  with  the  z  axis  as  the 
axis  of  symmetry,  the  stress-strain  relations  for  an  azimuthally-isotropic  medium  can  be  put  in  the 


convenient  form: 


T=  [C13  I  divu  +  c44(Vu  +  uV)]+  Ta 


Tania—  (cn  —  C13  —  2c44)(  /  — ezez)  +  ~Qy^)  (B-2) 

T  (C33  —  Cj3  2c44)-^— ezez 

,  ,  sL^u®.  -  ,  r>®uy  -  -  (®ux  &uy\,.  -  -  -  s' 

+  (c44  -  c66)  [2^eyey  +  2—exex  -  j  (exey  +  eyex)  . 

Here  ( ux,uy,uz )  are  the  cartesian  displacement  components  and  {cii,Ci3,C33,c44,c6r,}  are  the  five 
elastic  constants  of  the  medium.  The  unit  dyadic  is  /  and  {ex,ey,ez}  are  the  three  cartesian  unit 
vectors  along  the  respective  axes. 

We  now  take  the  divergence  of  T  in  (B.2).  Using  some  vector  and  tensor  identifies  (Ben- 
Menahem,  1987),  we  find  that  (B.l)  can  be  recast  in  the  form 


(ci3  +  c44)grad  div  u  -f  curl(ezL )  +  ezB  +  Vt  N  +  pu2ii  =  -pF  , 
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where 


L  =  -c66  A  -  c44  d]  {V72A}  , 

B  =  (C33  -  Ci3  -  C44)  +  C44  V2u2  , 

N  =  (cn  —  Cj3  —  C44)0  +  c44$j  20}  , 


A  =  e,  •  curl  u  =  %*■  -  , 


0  =  ^  +  9^  J 

v(  =  V  -  e2d2  =  exdx  +  evdy  ;  V2  =  V2  -  <92  .  (I 

and  V72  is  the  inverse  transverse  Laplace  operator,  to  be  given  explicitly  in  (B.33). 

We  apply  to  (B.3)  the  linear  operator  {e2  •  curl}.  This  yields  a  single  w'ave  equation  for  A 


c«vl\  +  ct<0] a  +  =  -e  (%  -  §7)  , 


where  (A,  Y.  Z)  are  the  cartesian  components  of  F. 

Next,  we  take  the  z-component  of  (B.3)  and  differentiate  it  with  respect  to  z.  The  result  is  a 
coupled  wave  equation  in  the  potentials  0  and  T 

A  <7 

[c44V2+O33O'2]r  +  (c13  +  C44)6>20  +  ^2r=  -p  —  .  (B.S) 

A  second  equation  is  obtained  from  (B.3)  by  taking  its  transverse  divergence 

[c4402  +  Cn^f]  ©  +  (ci3  +  C44 )A2 L  +  pu;20  =  —  p  ( +  ■— •  (B.9) 

We  shall  obtain  explicit  solutions  for  the  three  potentials.  Once  these  are  known,  the  components 
of  if  are  expressible  in  terms  of  these  functions: 

~  /dO  0\\  .  /tf©  l)\\ 

ur-V(  uy  —  Vt  u;-Jl,h.  (B.10) 


1 18 


Equation  (B.7)  can  readily  be  solved  in  closed-form.  Let  us  specify  the  body  force  distribution 
as  a  localized  force  of  magnitude  F(oj)  and  fixed  direction  e,  at  r  =  tq 


pF  =  F{u)e6{f-  r0)  , 


(B.ll) 


where  F  is  a  force  per  unit  mass,  but  F(u)  is  force  proper.  Assume 


A  =  —  Fe  •  curl(ezg)  , 


(B.12) 


where  g  is  an  unknown  function.  A  substitution  of  (B.12)  in  (B.7)  renders  at  once  the  equation  for 


cee  +  c44  d\  g  +  pu)2g  =  -6{f-  f0)  . 


(B.13) 


Its  exact  solution  is  the  scalar  Green’s  function, 


9(r,u)  = 


4tt  i?r/S//pa5 


(B.14) 


R  =  yf, A2  +  Vsh(z  -  zo)2 »  A2  =  (x  ~  xo)2  +  (y  -  yo )2 , 


2  ^66  1 2  ^ 

VSH  -  ~ —  ’  kSH  —  —  . 

C  44  0, 4 


(B.15) 


An  alternative  way  to  derive  (B.14)  is  instructive  at  this  point.  Assume  a  Fourier-integral 
representation  of  the  solution  of  (B.13) 

g(r^)=J^r  [  [  °)+Ky-yo)+M*-:o)]dad0d  (B16) 

Sn-ipJ-voJ  J  h(a,p,  7;w)  ' 

and  use  the  known  identity 

«5(r-r0)=^3  J°°  J  f  )]  dadpd-,  .  (B.17) 
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Inserting  these  into  (B.13),  yields 


/i(a,/?,7;w)  =  aH(a2  +  /r)  +  -5— a47‘ -  u  . 

Vsh 


(BAS) 


We  make  a  change  of  variables  from  (a,/?, 7)  to  spherical  coordinates  in  wave-number  space  (Fig¬ 


ure  2) 


UJ  -  -  L*J  ,  ~  ~  UJ 

a  =  —  sin  6  cos  <£  ;  /3  =  —  sin  9  sin  <j> ;  7  =  —  cos  9 

OOC- 


a2  +  /32  +  72  =  - 


3  — 

The  Jacobian  of  this  transformation  is  J  =  ^sin#  ,  and 


(B-19) 


/i(a,/3,7)=  ^(C2sh-C2), 


(B.20) 


where 


Csh  =  [«4  sin2  9  -F  a5  cos2  0]1/2  , 


(B-21) 


is  the  plane-wave  phase-velocity  of  the  SH  wave.  Substituting  (B.20)  and  (B.21)  into  (B.16),  the 
integral  becomes 


u>  y° 

8  x3p  Jo 


dC 

C2[C2  -  ClH{9)) 


[  sin  9d0e~'*zcosS  e-'fsine^co^+ysin 
Jc±  Jo 


(B.22) 


Since  7  =  -  (a2  +  /32),  where  (q2  +  /?2)  varies  from  zero  to  infinity,  7  can  take  on  pure 

imaginary  values.  In  isotropic  media  C±  is  the  usual  Weyl  contour  (Fig.  1).  However,  in  our 
anisotropic  medium,  there  are  two  branch  points  on  C*  at  r0  =  ±  tan  One  can  see  this 

as  one  puts  9  —  y  ±  it  into  the  expression  for  C§H(9)  via  Eq.  (B.21).  At  these  branch  points 
Csh  vanishes.  For  this  reason  the  integration  along  r  on  C*  must  stop  at  y  ±  tip.  As  775//  — *  1 
(isotropy)  tq  —*  ±00  and  the  original  Weyl  contour  is  re-established.  The  variation  of  7  is  then 


mapped  onto  a  corresponding  variation  of  6  along  either  of  the  modified  Weyl  contours  C+  or  C~ 
(Fig.  1),  depending  on  the  sign  of  z,  to  secure  finite  values  of  e~nz.  The  function  Cjw(0)  is  real 
for  all  complex  and  real  values  of  0  along  C±  ■ 

The  integral  over  4>  is  immediate  and  equals  2tJq  (^Asin#),  where  J0  is  the  Bessel  function  of 
the  first  kind  and  order  zero.  Then,  the  integral  over  C  is  evaluated  through  the  method  of  residues 
at  the  contributing  pole  C  =  Csh >  yielding 


/ 


e-t%zcoa6  jo  (£Asin  0)  XI  -,^_,cos.r  ,  „  A  .  - 

dcew^wr,  o]-cir 


-  z  COS  8  r  (  <*> 

0  Vi 


kCsh 


(B.23) 


The  net  result  of  both  integrations  over  <j>  and  C  is 

-tu  f  dOsin  0  (  u  -\  -t-*L-zcoae 

4x p  Jc±  C%h(9)  \Csh  J 

In  the  isotropic  limit  Csh  — *•  b.  Then  the  classical  Sommerf eld-  Weyl  integral 


(B.24) 


e-tkR 


-ikbR 


=  J  dOsinOJo  ^Asinflj  e 


cos  9 


(B.25) 


secures  the  reduction  of  g  to  its  isotropic  limit.  We  effect  another  change  of  variable 

*CSH  (f) 


.  .  \  2  /  n'-'sri  vov  .  .  sin  Odd 

smu  =  sin0— — ttt- ;  gsH  cos  u  =  cos  9— — -f-  ;  sin  u  du  =  — 


Csh{9) 


1  A  ■;  ;  sin  Udu=^-—-  Y  (B.26) 
Csh{0)  VSH  [Cs/H#)]3 

where  Cjw  (y)  =  04  and  ksH  =  ^7— •  One  can  easily  show  that  this  transformation  restores  the 
modified  Weyl  contour  into  the  standard  Weyl  contour  C±  (Fig.  1).  Therefore,  using  (B.25) 


-iLorjsH 


4xpa^2  Jc* 


J  sin  u  du  Jo(ksH^  sin  u)e 


l  e~ikSHR 


pgsHa  5  4  xR 


(B.27) 


in  agreement  with  (B.14). 

The  ensuing  displacements  are  obtained  from  (B.14)  via  (B.12)  and  (B.10) 

F(u) 


USH  = 


4  XTisHPas 


curl(e2x)  , 


(B.28) 
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X  =  curl(ezr )  •  e  ,  r  =  Vt2S ,  5(A,z)  = 


;  13.29) 


Next  we  use  the  tensor  identities  (Ben-Menahem,  1987) 


curl(G  ■ F )  =  (cur/  G)  •  F 


curl[ezcurl(ezr )]  =  /  -eze* - 5  , 

V,V,  1  [ 

— 2  $  =  eAe&S  +  (i^et  -  i£.eA)-£j  J  &S(A,z)dA  , 

i^5(A,2)  =  J^J&S(A,z)dA, 


( B.30 ) 


(B.32) 


(B.33) 


valid  for  a  fixed  vector  F,  a  variable  dyadic  G  and  remembering  that  S  is  independent  of  the 
azimuth  angle  d>.  Here  (ea,e^,e*)  are  the  unit- vector  triad  of  cylindrical  polar  coordinates.  We 


note  that  the  integral  in  (B.32)  renders 


J  AS(A,z)dA  =  -- 


_2 _ e~'ksHfl 


(B.34) 


Using  (B.32),  one  may  rewrite  (B.28)  as 


v-SH  = 


4irprisHas 


7  „  „  VtV,  A  ]  . 

I  -ezez  -  —x-  )S\-e 


=  F(u>)  Gsh  e  , 


(B.35) 


where 


7  x  1  (~  *  .  VtVt\e-fc™* 

Gstf  (ur0;u;  =  -  /  -e*e* =r~  ^ — 

^VSHpa*  \  V;  I  R 


(B.36) 


is  the  SH  part  of  the  total  Green’s  tensor.  Note  that  since  A  =  rsin0,  z  -  z0  =  r  cos  6  [0  is  the 


spherical  collatitude  angle  measured  from  the  z  axis],  we  have 


Ce6  2 


R  -  r  sim  6  -f  — —  cosJ  0 


(B.37) 
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Also  note  that  the  frequency  enters  into  the  displacements  only  via  the  source  function  F(u> )  and 

_i  —  R 

the  phase  e  . 
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APPENDIX  C.  LIST  OF  SYMBOLS 


a 


ai  =  cu/p 


P-wave  velocity  in  isotropic  media 


a2  =  C33/P 


a3  =  (C13  +  c44)/p  Square  of  velocities  in  Azimuthally-Isotropic  Media  (AIM). 
a4  =  C66/P 


^5  =  C44/p 


Ai  11  A 


& 


cijM 


c 


P? 


(Cm)2 


Cp 


C'sv 

Csh 


Plane-wave  response  tensor  (2.10). 

S-wave  velocity  in  isotropic  media  (=  y^. 

Tensor  of  elastic  coefficients. 

Matrix  of  elastic  coefficients  in  Voigt  abbreviated  notation. 

Weyl  contours. 

Phase- velocity  of  type-m  waves  (m  =  1,2,3)  in  general  anisotropic  media. 
Eigenvalues  of  Aij. 

Quasi-longitudinal  phase  velocity  in  AIM,  (m  =  1). 

Quasi-transverse  phase  velocity  in  AIM,  (to  =  2). 

SH-wave  phase  velocity  in  AIM,  (m  =  3). 

Unit  vector  in  the  direction  of  a  point- force  source. 

Unit  vectors  in  cartesian  coordinates. 
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^ A  >  ^ z 

?  € 8 't  €(f> 

F(t) 

F(u>) 

F 

9 

9ij 

G  ij ,  G 

Go 

Guo 

Gp 

Gsv 

Gsh 

H 

«-*♦ 

/;/ 

-Wl 

kp 

ksv 

ksn 


Unit  vectors  in  cylindrical  coordinates. 

Unit  vectors  in  spherical  coordinates. 

Time-function  of  applied  point-force  source  per  unit-volume. 
Fourier-transform  of  F(t). 

Applied  force. 

Scalar  SH  Green’s  function  for  AIM. 

Tensor  of  metric  coefficients. 

Spectral  Green’s  tensor  in  general  anisotropic  media. 
Singularity  tensor  for  AIM. 

Green’s  tensor  for  homogeneous-isotropic  elastic  solid. 

Quasi  P-wave  Green’s  tensor  for  AIM  at  far- field. 

Quasi  SV-wave  Green’s  tensor  for  AIM  at  far-field. 

SII  Green’s  tensor  for  AIM  at  far-field. 

Heaviside  unit-step  function. 

Unit  dyadic;  unit  matrix. 

Bessel  functions. 

u>/a\/2 . 

/  t/2 

“7  <  • 

/  1/2 
u>/aA'  . 
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k 

k 

m 

Mij 

h 

N(s) 

P-wave 

r(x,y,z) 

fo 


Wave-number  vector. 

Horizontal  wave-number. 

Index  of  wave- type  in  general  anisotropic  medium  (also  used  in  (4.65)). 
Christoffel  tensor  (2.12). 

Spectral  moment  of  dipolar  source. 

Unit  vector  in  k  direction. 

Unit  normal  to  slowness  surface. 

Abbreviation  for  the  quasi-longitudinal  wave  (Section  4.2). 

Coor’.r  us  of  field-point. 

Coordinates  of  source-point. 

Coordinates  of  field-point  in  AIM  metric. 

Magnitude  of  field-point  position  vector. 


R  =  \r-  f0| 
R 

Rp 

R. 


si-si  ,s2,s3) 


hi  h,  h 


SV-wave 


Effective  source-receiver  distance  in  AIM  for  SH-waves. 

Effective  source-receiver  distance  in  AIM  for  quasi  P-waves  (5.60). 
Effective  source-receiver  distance  in  AIM  for  quasi  SV-waves  (5.79). 
Slowness  vector  and  its  components  ( s  =  \J s\  4-  s^)- 
Unit  vectors  in  the  direction  of  the  (si,S2,S3)  respectively. 
Abbreviation  for  the  quasi-transverse  wave  (Section  4.2). 
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t 


Time 


Displacement  polarization  vector  of  type-m  waves. 
u  Displacement  vector. 

w  Group- velocity  vector  in  general  anisotropic  media. 

Wp(s)  Magnitude  of  group-velocity  of  quasi  P-waves  in  AIM. 

ws(s)  Magnitude  of  group-velocity  of  quasi  SV-waves  in  AIM. 

a,/3, 7  Direction  cosines  of  wave-number  vector. 

dm  Eigenvalues  of  Mi j. 

T  Auxiliary  potential  in  AIM. 

Tm  Slowness  surface  for  type-m  waves  (m  =  1,2,3). 

rp  Slowness  surface  in  AIM  for  quasi  P-waves. 

T,  Slowness  surface  in  AIM  for  quasi  SV-waves. 

6ij  Kronecker’s  delta. 

fi(r)  Dirac’s  delta-function  in  3D. 

d*  =  h  i  =  is 

(A,<t>,z)  Cylindrical  coordinates. 

(  Parameter  of  azimuthal-isotropy. 

nsn 


Coordinate  scaling  factor  for  SII  motion  in  AIM  (B.15). 


Coordinate  scaling  factor  for  quasi  P-motion  in  AIM  (5.44). 


71p 


Vs 


0 

(M) 

M) 

ki  ; 


P 

n 

da 

da 


a 


p 


^3 

* 

^3 


Coordinate  scaling  factor  for  quasi  SV-motion  in  AIM  (5.64). 
Auxiliary  potential  in  AIM. 

Spherical  polar  angles. 

Spherical  polar  angles  in  wave-number  space. 

Principal  curvatures  of  slowness  surface. 

Gaussian  curvature  of  type-m  slowness  surface. 

Lame  parameters  in  isotropic  elastic  solid. 

Uniform  mass  density. 

Solid  angle. 

Scalar  surface  element  in  Tm. 

Surface-element  vector  in  rm. 

Generalized  collatitudinal  angle  in  AIM  for  quasi  P-waves. 
Generalized  collatitudinal  angle  in  AIM  for  quasi  SV-waves. 
Fundamental  potential  in  AIM. 

Quasi  P-wave  potential  in  AIM. 

Quasi  SV-wave  potential  in  AIM. 

Angular  frequency. 


Transverse  gradient. 


Integration  operator. 
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LEGEND  OF  FIGURES 


Figure  1.  Path  of  integration  in  the  complex  6  plane  of  the  Sommerfeld- Weyl  integral  for  anisotropic 
media. 

Figure  2.  Cartesian  (a,/?, 7)  vs.  spherical  ( 0, 4>)  coordinates  in  wave-number  space.  On  the  same 
axes  we  also  show  cartesian  (x,y,z)  vs.  spherical  (r,0,d>)  in  real  space. 

Figure  3.  Relative  error  of  phase  velocities  for  four  earth  materials  as  a  function  of  the  collatitudinal 
(zenithal)  angle  (see  Table  2).  Top:  Quasi  P  velocity.  Bottom:  Quasi  SV  velocity. 

Figure  4.  Geometry  of  the  slowness  surface  in  a  generic  case  for  which  rjp  >  1. 

Figure  5.  Geometry  of  the  wave  surface  in  a  generic  case  for  which  r)p  >  1. 

Figure  6.  Exact  and  approximate  slowness  surfaces  for  four  earth  models  (see  Table  2). 

Figure  7.  Top:  Vertical  radiation  patterns  in  the  far-field  for  the  vertical  quasi-longitudinal  dis¬ 
placement  due  a  vertical  point  force  for  two  types  of  anisotropic  media.  Broken  lines  indicate 
isotropic  limit.  Force  is  directed  from  origin  toward  0°  along  the  positive  symmetry  z-axis  and  6 
is  measured  in  vertical  plane  positively  counterclockwise.  Bottom:  Same  as  above  for  the  vertical 
quasi -transverse  displacement  component. 

Figure  8.  Top:  Vertical  radiation  patterns  in  the  far-field  for  the  collatitudinal  quasi-longitudinal 
displacement  due  a  horizontal  point  force  in  the  x  direction.  Displacement  field  vanishes  in  the 
isotropic  limit.  Bottom:  Vertical  radiation  patterns  for  the  radial  quasi-transverse  displacement 
due  a  horizontal  point  force  in  the  x  direction.  Displacement  field  vanishes  in  the  isotropic  limit. 

Figure  9.  Top:  Vertical  radiation  patterns  in  the  far-field  for  the  vertical  quasi-longitudinal  dis¬ 
placement  due  to  an  explosion  point  source.  Broken  lines  indicate  isotropic  limit.  Bottom: 
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Vertical  radiation  patterns  for  the  vertical  quasi-transverse  displacement  due  to  an  explosion 
point  source.  Displacement  field  vanishes  in  the  isotropic  limit. 

Figure  10.  Vertical  radiation  patterns  in  the  far-field  for  the  SV/P  ratio  of  vertical  quasi-transverse 
to  quasi-longitudinal  components  of  displacement  field  due  to  an  explosion  point  source.  SV/P 
ratio  vanishes  in  the  isotropic  limit. 

Figure  11.  Top:  Vertical  radiation  patterns  in  the  far-field  for  the  horizontal-radial  quasi-transverse 
displacement  due  to  an  explosion  point  source.  Displacement  field  vanishes  in  the  isotropic  limit. 
Bottom:  Vertical  radiation  patterns  for  the  collatitudinal  quasi-longitudinal  displacement  due 
to  an  explosion  point  source.  Displacement  field  vanishes  in  the  isotropic  limit. 

Figure  12.  Top:  Vertical  radiation  patterns  in  the  far-field  for  the  collatitudinal  quasi-transverse 
displacement  due  to  an  explosion  point  source.  Displacement  field  vanishes  in  the  isotropic 
limit.  Bottom:  Vertical  radiation  patterns  for  the  radial  quasi-transverse  displacement  due  to 
an  explosion  point  source.  Displacement  field  vanishes  in  the  isotropic  limit. 

Figure  13.  A  comparison  of  ugp  radiation  patterns  in  the  far-field  due  to  two  sources  in  azimuthally 
isotropic  Timber  Mountain  tuff.  Left:  Strike-slip  dislocation  in  the  plane.  Right:  Explosion. 

Figure  14.  A  comparison  of  uz,  radiation  patterns  in  the  far-field  due  to  two  sources  in  azimuthally 
isotropic  Mesaverde  clayshale.  Left:  Strike-slip  dislocation  in  the  xz  plane.  Right:  Explosion. 
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LEGEND  OF  TABLES 


Table  1.  Kernels  of  the  Hankel-Transform  representation  of  displacement  fields  of  various  sources 
in  unbounded  azimuthally  isotropic  media. 

Table  2.  Structural  constants  of  some  azimuthally  isotropic  earth  materials  reported  by  Thomsen 
(1986). 

Table  3.  Cylindrical  components  of  displacement  field  and  radiation  patterns  of  point  sources  in 
unbounded  azimuthally  isotropic  media. 

Table  4.  Spherical  components  of  displacement  field  and  radiation  patterns  of  point  sources  in 
unbounded  azimuthally  isotropic  media. 
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EXPLOSION 
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TABLE  2.  Structural  constants  of  some  azimuthally-isotropic  earth  materials. 
(Density  in  g  •  cm-3;  Cij  in  units  of  1010  dyn  •  cm-2). 


(S) 

(t) 

(*) 

(*♦) 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

Cll 

A 

pa  1 

A  +  2/i 

50.78 

51.77 

50.01 

10.20 

3.52 

11.81 

56.91 

C33 

C 

pa-i 

A  +  2fi 

36.85 

48.29 

45.05 

7.02 

2.02 

9.68 

54.72 

Cl3 

F 

- 

A 

21.49 

-2.68 

-8.60 

4.95 

2.03 

7.11 

36.99 

<"66 

N 

pa  a 

M 

14.87 

25.97 

26.59 

2.31 

0.42 

2.26 

9.71 

C  44 

L 

pa5 

P 

11.01 

24.50 

24.57 

1.37 

0.27 

1.70 

8.03 

Cl3  +  C44 

pa  3 

A  +  p 

32.50 

21.82 

15.97 

6.32 

2.30 

8.81 

45.02 

P 

- 

- 

P 

2.56 

2.69 

2.87 

2.00 

1.80 

2.25 

2.33 

V 

1.000 

1.162 

1.030 

1.040 

1.298 

1.247 

1.153 

1.100 

V, 

1.000 

1.034 

0.891 

0.838 

0.741 

0.829 

0.921 

0.779 

*?P 

1.000 

1.163 

1.111 

1.202 

1.315 

1.362 

1.122 

1.072 

e  (percent) 

0 

18.9 

3.6 

5.5 

22.5 

37.1 

11.0 

2.0 

($)  Crystalographer’s  notation 
(f)  Stoneley’s  notation 
(*)  Musgrave-Buch wald’s  notation 
(**)  Isotropic  limit 


(1)  Mesaverde  clayshale  (Thomsen,  1986) 

(2)  Mesaverde  sandstone  (Lin,  1985) 

(3)  Mesaverde  sandstone  (Lin,  1985) 

(4)  ~og  Creek  shale  (Robertson  &  Corrigan,  1983) 

(5)  Eocene  Wills  Point  shale  (Robertson  &  Corrigan,  1983) 

(6)  Pierre  shale  (White  et  al.,  1983) 

(7)  Timber  Mtn.  tuff  (Shock  et  al.,  1974) 
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Stationary  phase  approximation.  (Notation  and  convention  as  in  Table  3.) 
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Computation  of  complete  waveforms  in  general  anisotropic 
media  -  Results  from  an  explosion  source  in  anisotropic  medium 

B.  Mandal  and  M.  N.  Toksoz 

Earth  Resources  Laboratory,  Dept,  of  Earth  Atmospheric  &  Planetary  Sciences,  Massachusetts 
Institute  of  Technology,  4 2  Carleton  St,  ES4-400,  Cambridge,  MA-02142,  USA. 

Submitted  to  Geophysical  Journal,  September,  1989 


SUMMARY 

An  algorithm  of  generating  complete  waveforms  in  general  anisotropic  media  (21  constants)  is  pre¬ 
sented.  Complete  waveforms  have  been  computed  including  a  point  source  in  anisotropic  medium 
and  the  reflection  and  transmission  properties  of  individual  interfaces.  A  recursive  scheme  of  scat- 
terer  operators  and  the  numerical  wavenumber  integration  technique  are  used.  The  6x6  system 
matrix  A  of  each  layer  is  derived,  and  each  element  of  A  is  expressed  in  terms  of  coefficients  of 
two  horizontal  wavenumbers.  This  form  of  system  matrix  A  reduces  the  computation  time  and 
uses  the  wavenumber  as  an  inner  loop  (vectorizing  loop).  The  complex  upper  Hessenberg  Ma¬ 
trix  and  QR  algorithm  are  used  to  compute  simultaneously  the  eigenvalues  and  eigenvectors  of 
A.  The  upgoing  and  downgoing  eigenvalues  are  separated  according  to  radiation  condition  and 
the  properties  of  z-component  Poynting  vector.  Filon  integration  is  applied  to  evaluate  the  inte¬ 
grals  over  two  horizontal  wavenumbers.  Since  wave  propagation  in  anisotropic  medium  is  a  three 
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dimensional  problem,  it  requires  an  enormous  number  of  computations.  To  reduce  the  compu¬ 
tation  time,  the  wavenumber  loop  is  used  as  the  deepest  loop.  An  efficient  computer  code  has 
been  developed  for  vectorization  in  the  super  computer.  Both  source  and  receivers  can  be  placed 
at  arbitrary  depths.  Three  component  synthetic  seismograms  are  computed  from  an  explosion  in 
the  anisotropic  medium.  Examples  are  calculated  for  half-space  multilayered  anisotropic  media. 
Theoretical  seismograms  show  that  radiation  patterns  are  strongly  affected  by  the  anisotropy.  For 
example,  an  explosion  in  an  anisotropic  medium  with  horizontal  axes  of  symmetry  generates  P, 
SV,  and  Rayleigh  waves,  all  whose  amplitudes  vary  with  azimuths,  and  significant  transverse  (SH) 
waves  at  certain  azimuths.  These  effects  could  be  explained  by  the  non-spherical  behavior  of  source 
due  to  the  direction  dependence  of  elastic  moduli.  The  waveforms  are  further  complicated  by  the 
reflection  and  conversion  phenomena  of  waves  propagating  in  anisotropic  media. 

Key  words:  general  anisotropy,  complete  waveforms,  explosion  source,  wavenumber  integration, 
vertical  fractures 

INTRODUCTION 

This  paper  presents  the  computation  of  complete  wavefield  in  general  anisotropic  media  due  to  a 
point  source.  It  also  investigates  the  effects  of  anisotropy  on  the  synthetic  seismograms.  In  re¬ 
cent  years  increasing  numbers  of  observational  studies  in  seismology  have  required  interpretations 
in  terms  of  anisotropic  elastic  properties  since  a  small  degree  can  significantly  affect  the  imaging 
of  the  wave  field.  This  variation  may  be  caused  by  a  variety  of  mechanisms,  including  miner¬ 
alization/deposition  orientations,  tectonic  stress  or  preferentially  oriented  fracture  patterns.  An 
anisotropic  medium  has  some  distinct  effects  on  data  which  cannot  be  explained  by  an  isotropic 
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or  heterogeneous  medium.  The  existence  of  anisotropy  in  the  crust  has  been  widely  detected  in 
field  data,  as  well  as  in  laboratory  experiments  (e.g.  Robertson  and  Corrigan,  1983;  Crampin  et 
al.,  1984;  Gaiser  et  al.,  1984;  Babuska,  1981).  The  theoretical  developments  of  wave  motion  in 
anisotropic  media  for  last  50  years  were  reviewed  by  Crampin  (1981)  and  Pao  (1983).  Recently, 
computation  of  synthetic  seismograms  in  anisotropic  media  have  been  developed  by  Booth  and 
Crampin  (1983)  and  by  Fryer  and  Frazer  (1984)  which  further  advances  toward  complete  waveform 
computation  in  general  anisotropic  media. 

To  compute  complete  synthetic  waveforms,  we  use  the  numerical  wavenumber  integration 
method  (e.g.  Bouchon  and  Aki,  1977).  In  this  method  the  computation  can  be  divided  into 
three  steps.  First,  the  steady-state  wave  field  radiation  by  a  seismic  source  is  represented  as  a 
superposition  of  waves  propagating  with  discrete  phase  velocities.  In  this  way  a  single  source  is 
replaced  by  a  periodic  array  of  sources.  The  effect  of  secondary  sources  is  removed  by  introducing 
both  a  small  constant  imaginary  part  of  the  frequency  and  proper  wavenumber  step  size.  Then  each 
phase  velocity  component  of  the  source  wave  field  is  propagated  through  a  stack  of  layers  from  the 
reflection  and  transmission  properties  of  individual  interfaces  using  a  recursion  scheme  of  scatterer 
operators  and  scatterer  products  (e.g.  Saastamoinen,  1980;  Kennett,  1983,  Fryer  and  Frazer,  1984). 
The  resulting  contributions  are  summed  for  each  wavenumber  to  transform  from  the  wavenumber 
domain  to  the  space  domain.  Finally,  the  solution  obtained  is  Fourier  transformed  to  the  time 
domain,  and  the  unwanted  effects  of  source  periodicity  and  complex  frequency  are  removed.  This 
final  solution  yields  the  complete  wavefield  within  the  medium  for  the  particular  source  considered 
over  a  specified  time  length.  This  method  is  common  for  isotropic  and  anisotropic  media  except  for 
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the  computation  of  the  eigensolutions.  In  an  isotropic  medium  or  transversely  isotropic  medium 
with  a  vertical  axis  of  symmetry,  this  computation  of  eigensolutions  can  be  performed  analytically 
(e.g.  Mandal  and  Mitchell,  1986).  The  energy  propagating  from  source  to  receiver  is  restricted  to 
the  sagittal  plane.  However,  in  azimuthal  anisotropic  media,  the  energy  generally  propagates  out  of 
the  sagittal  plane  causing  the  transverse  component  of  the  wavenumber  vector  to  be  non-zero  (e.g. 
Auld,  1973).  This  phenomenon  forces  us  to  compute  three-dimensional  eigensolutions  to  achieve 
complete  wavefields  in  anisotropic  media.  In  general,  the  elastodynamic  solutions  are  computed 
from  the  Christoffel  equation  of  six  order  polynomial  (e.g.  Fedorov,  1968;  Musgrave,  1970;  Frayer 
and  Frazer,  1987).  This  procedure  needs  a  high  degree  of  numerical  precession.  To  avoid  these 
problems,  we  derived  the  system  matrix  A  in  a  spatial  form  and  calculated  the  elastodynamic  so¬ 
lutions  by  computing  the  eigenvalues  and  eigenvectors  of  the  system  matrix  A  using  QR  algorithm 
(e.g.  Golub,  1983).  This  method  is  very  stable  and  eigenvalues  and  its  eigenvectors  are  computed 
simultaneously.  We  then  substitute  these  eigensolutions  in  the  reflectivity  and  numerical  wavenum¬ 
ber  integration  algorithm  to  get  the  complete  wave  field.  To  reduce  the  computation  time,  we  used 
Filon  quadrature  (1928)  relation  for  numerical  integration  in  wavenumber  domain.  This  technique 
will  allow  us  to  use  larger  wavenumber  step  size  and  smaller  oscillation  of  the  kernel  of  the  integral. 

In  this  computation  we  choose  to  use  a  Cartesian  coordinate  system.  Although  a  cylindrical 
coordinate  is  the  right  choice  to  compute  a  synthetic  wavefield  for  modeling  elastic  waves,  there 
is  no  simple  solution  for  general  anisotropic  media.  The  problem  of  wave  propagation  in  general 
anisotropic  media  is  thus  much  more  conveniently  addressed  in  Cartesian  coordinates,  with  the 
wavenumber  to  space  transformation  accomplished  by  numerical  integration  and  the  frequency  to 


time  transformation  by  inverse  Fourier  transformation. 

Investigated  is  the  complexity  of  seismic  radiation  patterns  and  the  generation  of  transversely 
polarized  (SH)  waves  by  an  explosion  in  anisotropic  medium.  The  causes  of  the  azimuthal  de¬ 
pendence  of  radiation  patterns  of  P,  SV,  and  Rayleigh  waves  and  the  generation  of  SH  and  Love 
waves  by  underground  nuclear  explosions  have  been  studied  extensively  (Kisslinger  et  al.,  1961; 
Archambeau  et  al. ,  1970;  Toksoz  and  Kehrei,  1972;  Masse,  1981;  Wallace  et  al.,  1983;  Gupta  and 
Blandford,  1983;  Lynnes  and  Lay,  1988;  Johnson,  1988).  Various  mechanisms  including:  tectonic 
strain  energy  release  by  relaxation  or  triggering  of  an  earthquake;  dislocation  across  cracks;  spalla¬ 
tion  and  “slapdown”;  and  scattering  from  heterogeneities  have  been  proposed  and  found  to  explain 
some  of  the  data.  No  single  mechanism  has  been  identified  that  could  satisfactorily  explain  all  the 
data  (Masse,  1981;  Gupta  and  Blandford,  1983;  Johnson,  1988.) 

The  evidence  for  seismic  anisotropy  in  the  earth’s  crust  has  been  also  increasing.  Recent  ob¬ 
servational  studies  showed  that  the  crust  material  may  have  seismic  anisotropy  due  to  aligned 
fractures  randomly  distributed  (e.g.,  Stephen,  1981;  Crampin,  1984;  Thomsen,  1986;  Winterstein, 
1986).  Several  major  factors  contribute  to  seismic  anisotropy  including:  (1)  preferred  orientation 
of  the  minerals  due  to  deposition  or  metamorphism;  (2)  geometric  effects,  such  as  alternating  high 
and  low  velocity  thin  beds  (e.g.,  shales,  carbonates);  and  (3)  stress-induced  preferred  orientation 
of  micro-  and  macro-fractures  in  shallow  crust.  The  complete  seismogram  synthesis  in  transversely 
isotropic  media  with  a  vertical  axis  of  symmetry  (Mandal  and  Mitchell,  1986)  showed  that  the  most 
simple  kind  of  anisotropy  can  alter  the  waveform  significantly.  The  reflection  and  transmission  phe¬ 
nomena  for  anisotropic  media  are  very  much  different  than  isotropic  case  (e.g.  Rokhlin  et  al.,  1986, 
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Mandal,  1989).  In  addition,  the  some  radiation  patterns  are  altered  if  the  source  is  located  ui  an 
anisotropic  medium.  For  example,  an  explosion  source  in  anisotropic  full  space  could  generate  a 
significant  amount  of  shear  wave  energy  where  none  would  be  generated  in  an  isotropic  medium 
(Mandal,  1988;  Ben-Menahem  and  Toksoz,  1989).  These  studies  also  show  that  the  radiation 
pattern  of  P  waves  from  an  explosion  in  a  transversely  isotropic  medium  is  not  spherical. 

In  this  paper  we  investigate  the  effect  of  anisotropy  on  seismic  radiation  by  computing  complete 
synthetic  seismograms  from  an  explosion  source.  We  take  a  simple  half-space  or  a  layer  on  a  half¬ 
space  where  the  axis  of  symmetry  is  horizontal  and  velocities  vary  with  azimuth.  Such  a  medium  can 
be  derived  from  aligned  vertical  fractures.  We  show  that  anisotropy  could  have  significant  effects 
on  radiation  patterns  of  all  phases  producing  a  transverse  component  of  P-wave  and  generating  a 
strong  SH-type  transverse  wave. 

THEORETICAL  BACKGROUND 
I.  Wave  Equation  and  Plane  Wave  Solution 

We  consider  a  layered,  anisotropic  media  bounded  above  by  a  free  surface  and  below  by  a  half-space 
(Figure  1).  The  material  properties  within  each  layer  is  homogeneous.  Any  sublayer  is  characterized 
by  a  thickness  H,  a  density  p ,  and  21  independent  stiffness  coefficients  with  a  6  x  6  (symmetric) 
matrix.  We  take  z-axis  positive  downward  and  assume  density  and  elastic  constants  are  functions 
of  z  only.  The  wave  equation  without  body  force  will  be: 

pdfUi  =  Tijj  =  (CijkiUk'i)j 


(1) 


where  displacements:  u  =  (ux,  uy,  uz)T,  stress  tensor:  r,j  =  Cijki^kl  (Generalized  Hook’s  Law), 
and  strain  tensor:  eu  ■ 

We  assume  a  displacement-stress  vector  b  whose  variables  are  continuous  across  any  horizontal 
plane,  and  whose  trial  plane  wave  solutions  is 


where  U  and  T  are  new  coefficients  of  the  displacement-stress  vector.  They  are  expressed  as, 

U  =  {Ux,Uy,Uz)T  coefficients  of  component  of  displacement,  r  =  ^j(txz,  Tyz,rzz)T  and  T  = 
-^j{TXZ,  TyZ,  Tzz)t  vertical  components  of  stress  tensor  and  its  coefficients.  The  px  and  py  are 
horizontal  phase  slownesses  and  u>  is  the  angular  frequency. 

The  substitution  of  these  trail  solutions  (2)  on  the  wave  equation  (1)  and  after  some  algebraic 
mesh  (Appendix  A),  we  will  end  up  a  set  of  linear  differential  equations  with  simple  form 
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that  they  behave  same  at  the  boundary.  Takeuchi  and  Saito  (1972)  used  same  approach  to  solve 
the  differential  system  for  transversely  isotropic  media. 

The  differential  system  (3)  can  be  solved  subject  to  specified  boundary  conditions  to  compute  the 
displacement-stress  field  (b)  at  any  depth  provided  the  medium  parameters  vary  with  z-direction. 
The  classical  propagator  matrix  method  can  be  used  to  transmit  the  wavefield  from  one  point  to 
other.  The  techniques  of  propagator  matrix  method  are  well  defined  (e.g.  Gilbert  and  Backaus, 
1966).  For  vertical  stratified  medium,  using  a  propagator  matrix  approach,  the  response  at  depth 
z  may  be  written 

b(z)  =  eiuAWx~Zo)b(z0) 

=  P(z,z0)b(z0),  (5) 

where  P(z,z0)  =  is  the  stress-displacement  propagator  (also  called  matricant  and 

matrizant  of  layer  matrix  for  a  homogeneous  medium)  at  depth  z  from  a  reference  depth  z0.  The 
detail  properties  of  the  propagator  has  been  given  in  many  previous  articles  (e.g.  Gilbert  and 
Backaus,  1966,  Woodhouse,  1974). 

From  the  definition  of  the  propagator  (e.g.  Gilbert  and  Backaus,  1966)  we  may  write 

dzP(z,z0)  =  iu>A(z)P{z,z0), 

P(z,z0 )  =  I,  (6) 

P(z,z0)  = 

where  I  is  a  6  x  6  identity  matrix;  D  is  the  eigenvector  matrix  of  A  (i.e. ,  each  column  of  D  is  an 
eigenvector)  and  A  is  the  diagonal  matrix  of  eigenvalues.  We  can  write  the  relation  between  D,  A, 


and  A  as, 


D-'AD  =  A.  (7) 

Some  properties  of  eigenvectors  ( D )  are  identified  provided  they  are  normalized  correctly  with 
respect  to  vertical  energy  flux  (e.g.  Chadwick  and  Smith  1977;  Garmany  1983;  Fryer  and  Frazer 
1984). 

II.  Eigensolutions 

A  vertically  stratified  medium  may  be  defined  by  a  series  of  homogeneous  layers.  To  propagate 
wavefield  from  one  point  to  another,  we  need  to  evaluate  eigenvalues  and  eigenvectors  of  the  system 
matrix  (A)  for  each  layer.  There  are  six  eigenvalues  with  six  corresponding  eigenvectors.  In 
general  anisotropic  medium,  major  portions  of  the  computation  has  been  devoted  to  compute 
the  eigenvalues  and  eigenvectors  of  the  system  matrix.  There  is  no  simplified  solution  for  this 
in  anisotropic  medium.  Hence,  one  has  to  go  for  numerical  solution  and  the  present  computer 
revolution  allows  us  to  do  that  job  efficiently.  There  are  some  practical  analytical  solutions  for  the 
medium  with  a  higher  degree  of  symmetry  (Mandal  and  Mitchell,  1986;  Fryer  and  Frazer,  1987). 

In  general,  the  eigenvalue  problem  in  anisotropic  medium  may  be  expressed  as  in  terms  of  wave 
normal 

(Ci:kin:m  -  pv26ik)Pk  =  0  (8) 

and  in  terms  of  slowness  vector 

(CijkiPjPl  -  pStk)Pk  =  0  (9) 

where  Pk  is  the  polarization  unit  displacement  vector  (eigen  vectors),  n  is  the  wave  normal,  and  v 
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is  the  phase  velocity.  These  equation  (8)  or  (9)  have  non  trivial  solutions  'or  Pk  #  0  only  if  the 
determinant  of  the  system  equals  zero,  know  as  Christoffel’s  matrix: 

\Cijkinjni  -  pv26ik\  =  0  (10) 

\CijkiPjPi  -  pSik\  =  0  (11) 

These  equations  (10-11)  are  used  to  determine  three  unknown  phase  velocities  and  corresponding 
vectors  of  the  permissible  elastic  waves  in  the  given  wave  direction  n. 

The  nature  of  the  wave  phenomenon  in  anisotropic  media  is  significantly  different  from  that  in 
isotropic  media.  The  properties  of  wave  phenomena  in  anisotropic  medium  are  discussed  elsewhere 
(  e.g.  Fedorov,  1968;  Musgrave,  1970;  Auld,  1973;  Crampin,  1981).  The  Christoffel  equation  (11) 
gives  6  roots  in  each  medium,  three  of  which  are  up-going  ( U )  and  three  of  which  are  down-going 
( D )  corresponding  to  three  fundamental  waves.  They  are  a  quasilongitudinal  wave  ( qP )  and  two 
quasitransverse  waves  (Si-  faster  and  S2-  slower  transverse  wave).  In  general,  none  of  these  will  be 
purely  longitudinal  or  purely  transverse  waves.  Each  of  the  waves  has  different  phase  and  group 
or  energy  (ray)  velocities.  The  group  or  energy  velocity  is  greater  than  or  equal  to  the  phase 
velocity  and  its  direction  does  not  coincide  with  the  wave  normal  (e.g.  Fedorov,  1968;  Auld,  1973, 
Rokhlin,  et  al.,  1968).  The  determination  of  the  eigenvalues  and  eigenvectors  using  equations  (10) 
and  (11)  have  been  reported  by  many  authors  (e.g.  Fedorov,  1968;  Musgrave,  1970;  Frayer  and 
Frazer,  1987;  Rokhlin,  et  al.,  1986).  To  solve  the  six  order  polynomial  (11),  a  very  high  degree 
precision  is  necessary  using  an  analytical  method  especially  when  waves  are  in  an  evanescent  regime. 
The  Laguerre’s  numerical  method  (Press  et  al.,  1986)  may  be  used  to  compute  the  roots  of  the 
polynomial  even  in  normal  precision  (Mandal,  1988). 
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Here,  we  are  using  a  simultaneous  method  of  computation  of  eigenvalues  and  eigenvectors  of 
system  matrix  A.  The  algorithm  of  this  method  is  taken  from  EISPACK  routines  (e.g.  Smith,  et  al, 
1976;  Golub  and  Van  Loan,  1983).  The  algorithm  has  few  steps.  First,  it  balances  the  6x6  complex 
matrix  (A)  to  reduce  the  overall  norm  of  the  matrix.  Second,  it  establishes  the  balanced  matrix 
in  an  upper  Hessenberg  form  (e.g.  Wilkinson  and  Reinsch,  1971).  Finally,  the  QR  algorithm  (e.g. 
Golub  and  Van  Loan,  1983)  is  used  to  compute  simultaneously  the  eigenvalue  and  eigenvectors  of 
the  upper  Hessenberg  matrix.  The  QR  transformation  preserves  the  upper  Hessenberg  form  of  the 
original  matrix.  This  method  is  very  stable  and  computes  faster  (0(n2))  than  the  general  matrix 
method  (0(n3))  (e.g.  Press,  et  al.,  1986). 

The  up-going  and  down-going  wave  field  may  be  separated  by  the  radiation  condition  for  u  >  0 

Im(qD )  >  0  and  Im(qu)  <  0,  (12) 

where  qD  and  qu  are  the  vertical  slownesses  or  eigenvalues  for  down-going  and  up-going  waves 
respectively. 

The  above-mentioned  radiation  conditions  are  not  sufficient,  especially  when  there  is  no  imag¬ 
inary  component  of  the  eigenvalues  for  propagating  modes  through  a  nonattenuative  medium. 
Another  set  of  necessary  conditions  should  be  available  to  make  sure  that  the  propagating  modes 
carry  energy  into  the  other  solid  medium  and  that  the  evanescent  modes  decay  into  the  other  solid 
medium  as  a  function  of  depth.  The  Poynting  vector  approach  is  found  to  be  more  appropriate  to 
separate  the  up-  and  down-going  waves  along  with  the  radiation  conditions  (e.g.  Musgrave,  1970; 
Henneke,  1972).  This  vector  is  in  the  direction  of  the  gradient  of  the  slowness  surface,  not  along 
the  normal  to  the  wavefront.  For  separating  up-  and  down-going  vectors,  the  ^-component  of  the 
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wave  vector  is  sufficient.  This  component  is  given  by 

Pz  =  -\{v*xTyt  +  v*yTxl  +  v*zTzz)  (13) 

where  the  u*’s  are  conjugates  of  the  respective  particle  velocities.  The  criteria  of  the  evanescent  and 
nonevanescent  modes  are:  (1)  a  negative  imaginary  part  for  the  evanescent  mode;  (2)  a  positive  real 
part  for  the  nonevanescent  mode.  These  characteristics,  along  with  the  other  radiation  conditions 
(Eq.  12),  will  ensure  separate  correct  up-  and  down-going  wave  fields. 

III.  Reflection  and  Transmission  Response 

Our  model  is  composed  of  a  series  of  uniform  anisotropic  layers  with  a  free  surface  2  =  0  and  a 
halfspace  at  z  >  zi  (Figure  1).  We  calculate  the  propagator  forms  using  reflectivity  coefficients 
as  describe  by  Kennett  and  Kerry  (1979)  and  Kennett  (1983).  We  compute  the  response  at  the 
uppermost  interface  and  then  adding  the  effects  of  successive  deeper  layers. 

Ignoring  the  free  surface  effects,  the  relation  between  the  stress-displacement  vector  at  depth  0 
and  zi  from  (5)  and  (6) 

b(zL)  =  D(z£/+)W(zl+,0+)T>_1(0+)6(0),  (14) 

where  +  (plus)-sign  refers  to  the  expression  just  below  the  interface.  W  =  eM*z.-°)A  js  known  as 
"wave- propagator’  and  it  depend  ;  on  the  difference  between  the  current  depth  zi  (half-space)  and 
0  (free-surface).  In  terms  of  the  wavevector  in  the  half-space,  the  relation  looks  like 

v(zL+)  =  W(zL+,0+)v(0+),  (15) 
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In  terms  of  the  up-  and  down-going  wavevector,  the  above  equation  can  be  expressed  as 


vu(zl+) 

vd(zl+) 


Wn  Jy12  ®t/(0+) 


W2l  W22  /  t  »D(0+) 


where  we  partition  the  wave  propagator  into  3x3  submatrices.  We  can  apply  the  obvious  boundary 
conditions  for  the  wavevector  at  the  free  surface  (up-going  energy  will  equal  to  down-going  reflection 
energy)  and  at  the  half-space  (no  up-going  energy)  as 


t?i/(0+)  =  Rdvd(  0+),  vd(zl+)  =  Tdvd(0+) 


where  R  and  T  are  reflection  and  transmission  matrices.  If  we  consider  the  source  originating 
inside  the  layer  medium  then  from  the  equations  (16  and  17),  all  the  reflection  and  transmission 


coefficients  will  be 


0  I 


To  Ru 


Wn  Wl2  1  I  Rd  Tv 


W2l  W22 


I  0 


We  can  write  a  scatterer  or  reflection  matrix  of  the  equation  as 


Tv  Rd 


Ru  To 


Wn  ~Wu  W12 

WnWtf  W22  -  W2X Wn1  W12 


(e.g.  Kennett,  1974;  Kennett  and  Kerry,  1979).  It  is  useful  to  present  equation  (19)  in  the  symbolic 


R  =  p(W) 


The  term,  p,  is  known  as  a  scatterer  operator  which  maps  the  propagator  W  into  scatterers  R. 


The  properties  of  the  operator  are  discussed  in  an  article  by  Saastamoinen  ( 19S0).  This  concept  of 
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a  "scatterer  product”  can  be  introduced  for  studying  the  scattering  properties  of  a  stack  of  layers 
which  may  be  anisotropic  media  (Saastamoinen,  1980;  Fryer  and  Frazer,  1984). 

Due  to  continuity  of  stress  and  displacement  at  an  interface,  the  wave  propagator  W  at  depth 
z,  takes  a  form 

W(z,+,zt-)  =  D-'izi+Mzi-).  (21) 

The  interface  scatterer,  which  includes,  the  reflection  and  transmission  effects  is  from  equations 
(20)  and  (21) 

R(Zl+,z,-)  =  p[W(zt+,zt-)\.  (22) 

The  total  response  of  the  interface  and  the  uniform  medium  by  using  the  scatter  product  (*)  is 
R(^|  —  ,  Zt-1  +  )  =  &(•?»'+,  zi~  )  *  ^i-I  +  )•  (23) 

Thus,  the  overall  response  of  the  layered  structure  can  be  computed  from  the  above  relations  by 
adding  s  iccession  of  layer-interface  scatterers  from  the  top  of  the  layer  as  going  downward  or  from 
half-space  going  upward,  so  that  for  step  i  we  have  a  form 

R(-r,+,  0+)  =  R(-?»+,  Zi-i  +)  ★  R(z,_i +,  0+ ).  (24) 

In  this  <  pproach,  there  are  no  growing  exponentials  (Saastamoinen,  1980)  so  that  we  will  get  a 
stable  r*  -ult  even  at  high  frequences. 

IV.  Point  Source  in  Anisotropic  Medium 

Our  main  goal  is  to  compute  the  seismic  response  at  several  receivers  located  in  other  arbitrary 
positions  from  an  arbitrarily  placed  source.  The  basic  method  used  in  this  case  for  a  single  receiver 


is  described  in  detail  by  Kennett  (1983,  Ch.  7).  A  Cartesian  coordinate  system  (x,y,z)  is  used  to 
compute  the  complete  wavefield  in  anisotropic  media.  Here,  we  do  not  consider  any  approximation 
of  propagation  in  anisotropic  media.  To  account  for  all  kinds  of  anisotropic  effects,  we  must  integrate 
both  horizontal  wavenumbers  (kx,ky)  to  achieve  full  response  from  a  point  source.  Attenuation  is 
introduced  by  specifying  imaginary  parts  to  the  elastic  constants  (Crampin,  1981). 

The  inclusion  of  a  source  in  an  anisotropic  medium  is  discussed  by  Frayer  and  Frazer  (1984). 
Choosing  a  Cartesian  coordinate  is  very  helpful  to  introduce  a  point  source  in  anisotropic  media. 
The  original  homogeneous  differential  equation  (3)  is  no-longer  valid  for  the  presence  of  any  source 
and  the  equation  (3)  will  take  an  inhomogeneous  form 

db 

=  t'wAb  +  F  (25) 

az 

where  F  =  (0,0,0,  fx,  fy,  fz)  is  a  source  term  containing  the  body  forces.  With  the  presence 

of  a  point  source  in  the  medium,  we  will  have  a  discontinuity  across  the  plane  containing  the  source 
(depth  z,).  A  point  source  at  r  =  (0,0,  z,)  can  be  represented  by  a  point  force  h  =  (hx,hy,hz)T 
and  a  symmetric  moment  tensor  M,  with  a  body  force  equivalent  given  by  (Purridge  and  Knopoff. 
1964) 

/  =  h6(r  —  r3)  for  simple  point  force,  (26) 

/  =  —M  ■  V<Kr  "  rs)  for  double  couple  source,  (27) 


where  6  and  V  are  the  delta  function  and  gradient  operator  respectively.  If  5  is  stress-displacement 
discontinuity  at  depth  z3,  then  an  equivalent  wave  vector  discontinuity  is 


(2S) 
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This  vector  S(za)  is  obtained  from  equations  (26-27),  by  the  Fourier  transform  over  x  and  t/,  in 
the  form  (e.g.  Frayer  and  Frazer,  1984) 


S(z,)  =  — 


1_ 

ILJ 


(  \ 
0 


0 

hx 

h,  ii 


for  point  force, 


(29) 


S(Zs)  = 


0 

pxM  XX  +  PyMXy 

PrM xy  T  PyMyy 
Px  M  zx  +  PyM 


'  0  ' 


+  A(z,) 


A/*, 

My, 

A/„ 


for  double  couple  source. 


(30) 


For  example,  for  an  explosion  source,  point  forces  h  will  vanish  and  only  diagonal  components  of 
equal  magnitude  moment  tensors  (A/rj,  Myy,  Mzz)  will  contribute  to  the  source.  This  discontinuity 
wave  vector  is  used  to  compute  the  radiation  pattern  for  different  kinds  of  sources  in  an  anisotropic 
medium. 


For  example,  in  general  anisotropic  medium,  the  source  vector  for  explosion  source  will  be  (from 


Appendix  A), 


For  example  in  orthorhombic  medium  (9  elastic  constants),  source  vector  will  be 


V.  Response  at  the  Receivers 


(31) 


(32) 


Computation  of  a  seismic  wavefield  in  anisotropic  media  involves  a  very  large  number  of  compu¬ 
tations.  Without  a  supercomputer  it  is  almost  impossible  to  compute  synthetic  wavefields.  An 
algorithm  of  generating  complete  seismograms  for  fast  computation  on  a  super  computer  (CRAY) 
is  used.  In  this  method,  Kennett’s  derivation  of  response  (Chapter  7,  Kennett,  1983)  at  the  receiver 
is  modified  to  express  a  better  structured  algorithm  for  multiple  receivers  (e.g.  for  VSP  types  of 
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geometry).  This  formulation  allows  us  to  compute  all  the  matrices  in  a  single  pass  through  the 
layer  loop  and  helps  to  increase  the  speed  of  computation  in  vector  and  parallel  computers.  New 
relations  at  the  receiver  are  as  follows: 

for  receiver  depth  zr  <  za 

u(zr)  =  (. Mur  +  Mdrr(jR)  (Tl*yl  Tfr  (/  -  RU’vY  \rD^d{z,)  +  2u(z,j\  ,  (33) 

and  for  zr  >  za 

u(zR)  =  ( Mdr  +  MurR%l)  (l-RpRjpy1  Tg3  (i-r£r$)-1  \RfrZv{za)  +  ED(z,)]  ,  (34) 

where  RpL  is  expressed  in  the  algorithm 

nRL  _  rpaR  (  nsL  u sR\  ~ *  rpsR  .  nsR 
kD  —  1 D  \Md  ~  kD  )  lU  +  My 

These  notations  are  explained  in  Kennett  (1983).  These  modified  relations  (33-34)  are  also  used 
by  Mallick  and  Frazer  (1988). 

VI.  Numerical  Integration  and  Time  domain  Response 

Finally,  we  get  a  time  domain  response  due  to  a  point  source  at  the  receiver  by  the  relation 

1  f+°°  f+ °°  f+oo 

u{x,y,zR-,t)  =  /  /  /  u(kx,kv,zR,u)ei^x+k^-^dkxdkydu,  (35) 

[tvi )  y-oc  7~oo  ‘/-oc 

where  u(kx<ky,ZR,u),  the  displacement  (Green’s  function)  of  the  entire  medium  at  the  receiver,  is 
computed  for  each  kx,ky ,  and  u>  (from  relation  33-34)  and  for  the  algorithm  in  Figure  2.  Denoting 
by  ur  and  uj,  the  real  and  imaginary  parts  of  the  frequency,  the  impulse  response  u(x,y,ZR\  ()  can 
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then  be  written  in  terms  of  a  wavenumber  summation  form  (e.g.  Bouchon  and  Aki,  1977) 


r 


u(x,y,ZR;t)  = 


p— u>/t  r+oo  Ajr 

/  e,UJRtduR.-—~  m  u(^n  ,kym,zR\u)e' 

00  L*LV  l n=-Nm=-M 


e'kymV  ,  (36) 


where  Lx  and  Ly  are  the  secondary  source  intervals  along  the  x-  and  y-directions  and 


,  2tt  2tt 

kxn  =  7 -n  and  kym  =  — m  . 


The  A  and  M  are  wavenumber  loops  and  their  values  will  depend  on  the  convergence  of  the 
integral  kernel  u(kx,ky,  zr\uj).  The  imaginary  part  of  the  frequency  is  used  to  avoid  the  influence 
of  singularities  of  the  integration  kernel  (Bouchon,  1981).  We  later  remove  the  imaginary  part  of 
the  frequency  from  the  time-domain  solution. 

This  numerical  summation  (36)  is  only  applicable  to  the  ranges  ( R )  where  the  product  of  the 
kernel  and  the  exponential  function  is  well  represented  by  a  linear  function.  To  ensure  these  criteria, 
the  wavenumber  sampling  has  to  satisfy  the  following  inequalities  (e.g.  Bouchon  and  Aki,  1977; 
Mallick  and  Frazer,  1987), 


A kR  <  j,  R  <  ^  and  \J(L  -  Rf  +  z2  >  Vmaxt. 


Vmox  is  the  fastest  velocity  of  the  media.  For  larger  distance,  the  wavenumber  step  size  will  be 
smaller  and  hence  computation  time  will  increase.  However,  to  overcome  this  problem,  we  can 
solve  wavenumber  integration  by  using  Filon’s  quadrature  formula  (e.g.  Filon,  1928;  Frazer  and 
Gettrust,  1984).  In  this  method,  we  can  use  a  large  step  size  which  reduces  the  computational  time. 
Some  important  features  of  Filon’s  quadrature  formula  had  been  discussed  in  previous  articles  (e.g. 
Frazer  and  Gettrust,  1984;  Mallick  and  Frazer,  1988). 
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Since  we  use  the  discrete  wavenumber  summation  to  substitute  a  wavenumber  integration,  we 
need  to  compute  the  eigensolutions  and  reflection  and  transmission  matrices  of  each  medium  for  a 
large  number  of  wavenumbers  for  each  frequency.  We  use  wavenumber  (A:)  as  an  inner  loop  with 
the  main  intent  to  vectorize.  At  first,  we  compute  the  maximum  value  of  wavenumber  which  is 
sufficient  to  explain  all  types  of  waveform  within  a  given  time  window  for  each  frequency.  We  ' 

determine  the  total  number  of  wavenumber  steps  (JVfcmai)  for  the  proper  wavenumber  step  size. 

These  Ajtmax-loops  are  again  subdivided  into  sets  of  small  fc-loops  which  can  be  a  variable  according 
to  the  machine  memory  requirement.  This  code  can  be  executed  on  any  size  computer.  However, 
each  small  jfc-loop  vectorizes  well  on  a  super  computer  such  as  a  CRAY.  A  schematic  flow  chart  of 
this  algorithm  in  Figure  2  describes  the  flow  of  this  algorithm. 

In  summary,  the  depth-dependent  Green’s  function  is  computed  for  each  frequency  and  wavenum¬ 
ber  using  reflectivity  technique,  and  the  transformation  from  wavenumber  domain  to  space  domain 
is  performed  by  numerical  integration  which  requires  both  truncation  and  discretization  of  the 
wavenumber  interval.  Finally,  the  complete  response  in  a  time  domain  has  been  achieved  by  in¬ 
verse  Fourier  transformation  of  the  product  of  frequency  domain  Green’s  function  with  required 
source-time  function. 

We  used  a  simple  crustal  model  for  testing  this  algorithm  with  the  isotropic  model.  The  veloc¬ 
ities  of  the  top  layer  Vp  =  6.15  km/sec  and  Vs  =  3.55  km/sec  and  density  p  =  2.8  gm/cc.  The 
velocities  of  the  half-space  Vp  =  8.09  km/sec  and  Vs  =  4.67  km/sec  and  density  p  =  3.3  gm/cc. 

After  the  responses  are  computed  by  wavenumber  and  frequency  integration  with  0-2  Hz  band,  the 
synthetic  seismograms  are  generated  by  convolving  it  with  a  Ricker  wavelet-type  pulse  centered  at 
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1  Hz.  The  comparisons  of  synthetic  seismograms  at  a  distance  of  75  km  for  an  explosion  source 
depth  of  10  km  are  shown  in  Fig.  3.  The  synthetic  seismograms  on  the  left  are  representing  the 
vertical  (top)  and  radial  (bottom)  components  for  the  isotropic  computation  with  cylindrical  co¬ 
ordinates  (Mandal  and  Mitchell,  1986).  The  waveforms  in  the  right  represent  the  results  for  the 
same  medium  using  anisotropic  algorithm  of  this  article.  This  comparison  shows  that  the  present 
algorithm  of  general  anisotropic  medium  provides  results  which  are  consistent  with  other  methods 
for  a  model  with  isotropic  elastic  parameters. 

SYNTHETIC  SEISMOGRAMS  FROM  AN  EXPLOSION  SOURCE 
IN  ANISOTROPIC  MEDIUM 

To  understand  the  effects  of  anisotropy  on  seismic  waves  from  an  explosion  source,  we  use  two 
relatively  simple  test  models:  an  anisotropic  half-space  and  an  anisotropic  layer  over  an  isotropic 
half-space.  For  an  anisotropic  half-space,  we  investigate  the  effects  of  two  different  anisotropic 
degrees  of  the  anisotropy. 

I.  Half-Space 

The  anisotropic  half-space  is  our  simplest  test  model.  The  anisotropy  is  assumed  to  be  due  to  a 
parallel  set  of  thin  vertical  fractures  filled  by  a  fluid.  The  aspect  ratio  of  the  fractures  is  0.001.  The 
degree  of  anisotropy  of  the  models  is  determined  by  the  fracture  densities.  The  two  models  we  use 
have  fracture  densities  of  5%  and  20%.  The  elastic  constants  are  computed  using  the  theoretical 
formulation  of  Hudson  (1981)  and  Crampin  (1984).  Properties  of  the  anisotropic  medium  with 
parallel  vertical  fractures  are  described  by  a  set  of  five  elastic  constants  with  a  horizontal  symmetry 
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axis.  Models  for  estimating  the  effects  of  cracks  have  been  developed  (Hudson,  1981;  Schoenberg 
and  Douma,  1988;  Douma,  1988).  We  use  Hudson’s  (1981)  formulation  for  small  aspect  ratio 
fractures.  We  also  assume  fractures  are  fluid-filled.  Anisotropic  velocities  are  affected  by  fluid 
properties  such  as  air,  water,  etc.  (Hudson,  1981;  Shearer  and  Chapman,  1989). 

The  five  elastic  constants  in  terms  of  velocities  km/ sec )  are: 

1.  5%  fracture  density: 

un  =  6.12  V12  =  U33  =  6.15  v\i  =  «i3  =  3.53 

V44  =  3.55  t>55  =  i>66  =  3.35 

2.  20%  fracture  density: 

Un  =  5.94  V22  =  U33  =  6.13  V12  =  t)i3  =  3.43 

V44  -  3.55  v55  =  vee  =  2.85 

The  density  (p)  is  2.8  gm/cc  in  both  cases. 

An  explosion  source  is  placed  at  a  depth  of  1  km  below  the  surface.  Let  us  assume  the  fracture 
planes  are  parallel  to  an  east-west  direction  as  shown  in  Fig.  4,  and  the  Y  axis  (North)  is  normal 
to  the  fracture  plane.  We  calculate  seismograms  for  five  different  azimuths  (NO0,  N22.5°E,  N45°E, 
N67.5°E  and  N90°E)  with  stations  at  distances  5  and  50  km  from  the  epicenter.  After  the  responses 
are  calculated  by  wavenumber  and  frequency  integration  in  0-10  Hz  frequency  band,  the  synthetic 
seismograms  are  generated  by  convolving  these  with  a  Ricker  wavelet-type  pulse  centered  at  1  Hz, 
2  Hz,  3  Hz,  and  4  Hz.  The  seismograms  are  shown  in  Figures  5-6  for  5%  and  Figures  7-8  for 
20%  fracture  densities,  respectively,  for  the  five  azimuths  (<£)  in  the  0°  to  90°  range.  Prominent 
transverse  component  motion  is  observed  at  azimuths  of  22.45  and  67°.  At  5  km  offsets,  the  arrival 
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times  between  longitudinal  and  transverse  waves  are  small.  In  the  high  frequency  seismograms 
(4  Hz),  wave  types  can  be  differentiated  by  comparing  transverse  components  with  vertical  or  radial 
components.  In  this  and  subsequent  figures,  the  seismograms  in  the  direction  of  the  symmetry  axis 
(<p  =  90°)  '■orrespond  exactly  to  that  of  an  isotropic  medium  with  appropriate  velocities.  As 
the  degree  of  anisotropy  increases,  the  intensity  of  the  transverse  component  also  increases.  In 
an  anisotropic  medium,  the  three  wave  types  are:  qP  (quasi-longitudinal),  qSP  (transverse  wave 
polarized  within  the  symmetry  plane),  and  qSR  (transverse  wave  polarized  perpendicular  to  the 
symmetry  plane). 

The  seismograms  at  a  distance  of  50  km  from  the  source  (Figure  6)  show  that  qP  waves  show 
a  variation  of  travel  times  (velocities)  and  amplitudes  as  a  function  of  azimuth  on  both  vertical 
and  radial  components.  There  is  a  small  transverse  component  of  P  wave  at  <f>  =  22.5°,  45°,  and 
67.5°  azimuths  for  20%  fracture  density  (Figure  8).  The  large  pulse  in  the  vertical  component  is 
a  combination  of  shear  (qSP)  and  Rayleigh  waves.  The  velocity  variation  with  azimuth  is  clearly 
observable  as  is  the  amplitude  variation.  The  transverse  component  seismograms  clearly  show  the 
large  transverse  shear  wave  (qSR)  at  <f>  =  67.5°for  20%  fracture  density.  Also,  the  shear-wave¬ 
splitting  diagnostic  of  an  anisotropic  medium  is  observable. 

The  frequency  dependence  of  compressional,  shear,  and  Rayleigh  wave  amplitudes  may  be 
explained  by  the  depth  of  the  source  and  the  free  surface  effects.  In  seismograms  at  90°  azimuth 
(i.e.,  the  isotropic  case),  we  see  that  the  Rayleigh  wave  amplitudes  decrease  exponentially  with 
frequency.  At  other  azimuths  (22.5°,  45°,  67.5°),  the  qSR  and  Rayleigh  wave  amplitudes  change 
rapidly.  These  changes  are  less  for  5%  fractured  medium  (Figure  6).  In  anisotropic  media,  reflection 
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of  waves  from  an  interface  is  more  complex  than  for  the  case  of  isotropic  m^dia.  For  example,  there 
is  a  cusp  on  the  shear  wave  velocity  surface  (qSP)  for  fluid-filled  fractures  (Crampin,  1984;  Shearer 
and  Chapman,  1989).  Other  complexities  of  reflection  and  transmission  properties  in  an  anisotropic 
medium  have  been  described  elsewhere  (e.g.,  Rokhlin  et  al.  ,  1986;  Mandal,  B.,  1989). 

The  generation  of  transverse  components  by  an  explosion  in  an  anisotropic  half-space  is  due  to 
azimuthal  variation  of  elastic  constants  which  results  in  non-spherical  displacements  at  the  source. 
Although  we  construct  the  explosion  source  using  three  equal  orthogonal  moment  tensors,  the  final 
source  vector  is  a  product  of  the  moment  tensor  matrix  and  the  system  matrix  whose  components 
contain  medium  elastic  constants  ( Cijki ),  density,  and  phase  slownesses. 

II.  Layer  Over  Half-Space 

This  model  consists  of  a  2  km  thick  anisotropic  layer  over  an  isotropic  half-space.  The  source 
is  an  explosion  at  1  km  depth  in  the  anisotropic  layer,  as  shown  in  Figure  9.  The  anisotropic 
layer  contains  aligned,  vertical,  fluid-filled  fractures  with  an  aspect  ratio  of  0.001  and  fracture 
density  20%,  the  same  as  in  the  previous  model.  The  five  elastic  constants  in  terms  of  velocities 
=  are:  wn  =  5.35,  V22  —  w33  =  5.47,  vjj  =  V13  =  3.4,  V44  =  3.00,  v 55  =  =  2.42. 

The  density  (p)  is  2.7  gm/cc.  The  isotropic  half-space  has  the  P-  and  S-wave  velocities  of:  Vp  =  6.2 
and  Vs  =  3.5  km/sec.  The  density  is  2.8  gm/cc. 

Synthetic  seismograms  are  calculated  for  the  same  source  and  receivers  arrangements,  as  in 
the  previous  half-space  example.  Figures  10  and  11  show  the  synthetic  seismograms  at  distances 
5  km  and  50  km,  respectively.  The  transverse  motion  due  to  anisotropy  is  very  prominent.  The 
seismograms  are  more  complex  due  to  the  presence  of  multiple  reflections  and  normal  modes. 
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The  phenomena  for  anisotropic  media  have  some  very  different  features.  The  behavior  of  polar¬ 
izations  of  the  three  fundamental  waves  (qP,  qSP  and  qSR)  in  vertical  fluid-filled  fracture  medium 
has  been  studied  by  Shearer  and  Chapman  (1989).  These  computations  of  complete  wave  synthesis 
show  us  many  features  of  anisotropic  wave  propagation. 

CONCLUSIONS 

The  reflection  matrix  and  the  numerical  wavenumber  integration  are  used  to  compute  complete 
wavefields  in  stratified  anisotropic  media  due  to  a  point  source.  The  point  source  may  be  an 
explosion,  point  force,  or  double  couple  source.  The  source  and  receivers  can  be  placed  at  any 
arbitrary  place  in  the  medium. 

An  explosion  in  an  anisotropic  medium  has  a  non-uniform  radiation  pattern  and  generates  sig¬ 
nificant  transverse  motion  at  certain  azimuths  because:  (1)  A  spherical  explosion  in  an  anisotropic 
medium  becomes  non-spherical  due  to  the  directional  dependence  of  elastic  compliance.  (2)  The 
reflection  (especially  at  free  surface)  and  conversion  of  seismic  waves  propagating  in  an  anisotropic 
medium  further  complicates  the  seismograms  and  their  frequency  dependence.  These  examples 
show  that  anisotropy  may  be  a  significant  factor  in  contributing  to  the  generation  of  transverse 
waves  by  explosions.  From  these  synthetic  computations,  it  is  clear  that  anisotropy  takes  a  major 
roll  in  wave  propagation. 

The  complete  wavefield  in  anisotropic  media  is  now  possible  with  the  help  of  a  super  computer 
and  an  efficient  computational  algorithm.  We  hope  that  the  high  speed  super  computer  and  efficient 
code  will  give  us  the  opportunity  to  study  wave  propagation  in  general  anisotropic  media. 


175 


Acknowledgments 


This  work  is  partially  supported  by  the  Reservoii  Delineation  -Vertical  Seismic  Profiling  Consor¬ 
tium  at  M.I.T.  Partial  support  was  provided  by  the  Defence  Advanced  Research  Projects  Agency 
through  contract  F19628-87-K-0054  administered  by  the  Air  Force  Geophysical  Laboratory.  Some 
of  the  computations  have  been  done  on  the  CRAY  at  the  Pittsburgh  Supercomputing  Center  under 
NSF  Grant  EAR-8804355. 

References 

Alterman,  Z.,  Jarosch,  H.,  Pekeris,  C.  L.,  1959.  Oscillations  of  the  earth.  Proc.  Roy.  Soc.  Str.. 
A 252,  80-95. 

Archambeau,  C.G.  and  C.  Sammis,  1970.  Seismic  radiation  from  explosions  in  prestressed  media 
and  the  measurement  of  tectonic  stress  in  the  Earth,  Rev.  (Jcophys.,  6,  241  288. 

Auld,  B.A.,  1973.  Acoustic  Fields  and  Wave  in  Solids,  Vol.  1,  Wiley,  New  York. 

Babuska,  V.,  1981.  Anisotropy  of  vp  and  v,  in  rock-forming  minerals,  J.  Ccophys.,  50,  1  G. 
Ben-Menahem,  A.  and  M.N.  Toksoz,  1989.  Radiation  patterns  from  explosions  in  anisotropic  media, 
submitted  to  Science. 

Booth,  D.C.  and  Crampin,  S.,  1983.  The  anisotropic  reflectivity  technique:  theory,  Cnophys.  J.  R. 
astr ,  Soc.,  72,  755-766. 

Bouchon,  M.,  1981.  A  simple  method  to  calculate  Green's  functions  for  elastic  layered  media.  Hull 
Seism.  Soc.  Am.,  73,  959-971. 

Bouchon,  M.  and  K.  Aki,  1977.  Discrete  wavenumber  representation  of  seismic  source  wave  field-.. 


Bull.  Seism.  Soc.  Am.,  67,  259-277. 

Burridge,  R.  and  KnopofF,  L.,  1964.  Body  force  equivalents  for  seismic  dislocations,  Bull,  seism. 
Soc.  Am.,  54,  1875-1888. 

Chadwick,  1’.  and  Smith,  G.  D.,  19/7.  Foundations  in  the  theory  of  surface  waves  in  anisotropic 
elastic  materials,  in  Advances  in  Applied  Mechanics,  17,  ed.  Yih,  C.  S.,  Academic  Press,  New 
York. 

Crampin,  S.,  1984.  Effective  elastic  constants  for  wave  propagation  trough  cracked  solids,  Geophys. 
J.  R.  astr.  Soc.,  78,  691-710. 

Crampin,  S.,  1981.  A  review  of  wave  motion  in  anisotropic  and  cracked  elastic- media,  Wave  Motion, 
8,  343  391. 

Crampin,  S.,  E.M.  Chesnokov  and  R.G.  Hipkin,  1984.  Seismic  anisotropy— the  state  of  the  art:  II, 
Geophys.  J.  R.  astr.  Soc.,  76,  1-16. 

Douma,  J.,  1988.  The  effect  of  the  aspect  ratio  on  crack-induced  anisotropy,  Geophy.  Prospecting, 
36,  614-632. 

fedorov,  f  .1.,  1968.  Theory  of  Elastic  Waves  in  Crystals,  Plenum  Press,  New  York. 

1  ilon,  L.N .G.,  1928.  On  a  quadrature  method  for  trigonometric  integrals,  Proc.  Roy.  Soc.  Edinburgh, 
49,  38  47. 

1  razer,  L.N.  and  J.F.  Gettrust,  1984.  On  a  generalization  of  Filon’s  method  and  the  computation 
of  oscillatory  integrals  of  seismology,  Geophys.  J.  R.  astr.  Soc.,  76,  461  481. 

fryer,  G.J.  and  L.N.  Frazer,  1984.  Seismic  waves  in  stratified  anisotropic  media,  Geophys.  J.  R. 
astr.  Soc.,  78,  691  710. 


177 


Fryer,  G.J.  and  L.N.  Frazer,  1987.  Seismic  waves  in  stratified  anisotropic  media  -  II:  Elastodynamic 
eigensolutions  for  some  anisotropic  systems,  Geophys.  J.  R.  astr.  Soc.,  81,  73-101. 

Gaiser,  J.E.,  R.W.  Ward,  and  J.P.  DiSiena,  1984.  Three-component  vertical  seismic  profiles:  Velocity 
anisotropy  estimates  from  P  wave  particle  motion,  Handbook  of  Geophysical  Exploration,  section 

I,  Vol.  14B,  189-204. 

Garmany,  J.,  1983.  Some  properties  of  elastodynamic  eigensolutions  in  stratified  media,  Geophys. 

J.  R.  astr.  Soc.,  75,  565-569. 

Gilbert,  F.  and  G.E.  Backaus,  1966.  Propagator  matrices  in  elastic  wave  and  vibration  problems, 
Geophysics,  81,  326-332. 

Golub,  G.  H.  and  Van  Loan,  C.  F.,  1983.  Matrix  Computations ,  The  Johns  Hopkins  University 
Press,  Baltimore. 

Gupta,  I.N.  and  R.R.  Biandford,  1983.  A  mechanism  for  geneiation  of  short-period  transverse 
motion  from  explosions,  Bull.  Seism.  Soc.  Am.,  73,  571-591. 

Johnson,  L.R.,  1988.  Source  characteristics  of  two  underground  nuclear  explosions,  Geophys.  J.,  95 
15-30. 

Henneke,  E.G.,  1972.  Reflection-refraction  of  stress  wave  at  a  plane  boundary  between  anisotropic 
media,  J.  Acoust.  Soc.  Am.,  51,  210-217. 

Hudson,  J.  A.,  1981.  Wave  speeds  and  attenuation  of  elastic  waves  in  material  containing  cracks, 
Geophys.  J.  Roy.  astr.  Soc.,  64,  133-150. 

Kausel,  E.,  1986.  Wave  propagation  in  Anisotropic  layered  media,  Inter.  J.  Num.  Methods  in  Eng., 
23,  1567-1578. 


178 


Kennett,  B.  L.  N.,  1974.  Reflections,  rays,  and  reverberations,  Bull,  seism.  Soc.  Am.,  64,  1685-1696. 

Kennett,  B.L.N.,  1983.  Seismic  Wave  Propagation  in  Stratified  Media,  Cambridge  University  Press. 

Kennett,  B.  L.  N.  and  Kerry,  N.  J.,  1979.  Seismic  waves  in  a  stratified  half  space,  Geophys.  J.  R. 
astr,  Soc..  57,  557-583. 

Kisslinger,  C.,  E.J.  Mateker,  Jr.,  and  T.V.  McEvilly,  1961.  SH  motion  from  explosions  in  soil,  J.  of 
Geophys.  Res.,  66,  3487-3496. 

Lynnes,  C.S.  and  T.  Lay,  1988.  Analysis  of  amplitude  and  travel-time  anomalies  for  short-period 
P-waves  from  NTS  explosions,  Geophys.  J.,  92,  431-443. 

Mallick,  S.  and  L.N.  Frazer,  1987.  Practical  aspects  of  reflectivity  modeling,  Geophysics,  52,  1355- 
1364. 

Mallick,  S.  and  L.N.  Frazer,  1988.  Rapid  computation  of  multioflset  vertical  seismic  profile  synthetic 
seismograms  for  layered  media,  Geophysics,  53,  479-491. 

Malvern,  L.F.,  1969.  Introduction  to  the  Mechanics  of  a  Continuous  Media,  Prentice-Hall,  Engle¬ 
wood  Cliffs. 

Mandal,  B.,  1989.  Reflection  and  transmission  properties  of  elastic  waves  on  a  plane  interface  for 
general  anisotropic  media,  annual  report,  Reservoir  Delineation  Consortium,  ERL,  MIT.  2-1,2- 
48. 

Mandal,  B.,  1988.  Computation  of  the  complete  wavefield  in  anisotropic  media,  annual  report. 
Reservoir  Delineation-  Vertical  Seismic  Profiling  Consortium.  ERL,  MIT.  2-1,2-31. 

Mandal,  B.  and  B.J.  Mitchell,  1986.  Complete  seismogram  synthesis  for  transversely  isotropic  me¬ 
dia,  J.  Geophys.,  59,  149-156. 


179 


Masse,  13. 1’.,  1981.  Review  of  seismic  source  models  for  underground  nuclear  explosions,  Bull.  Seism. 
Soc.  Am..  71,  12-19  12G8. 

Musgrave,  M.J.P.,  1970.  Crystal  Acoustics,  Holden-Day,  San  Fransisco. 

Pao.  Y,  1983.  Elastic  waves  in  solids.  Journal  of  Applied  Mechanics ,  50,  1 152-1  Hi  t. 

Press,  W.H.,  13. P.  Flannery,  S.A.  Teukolsky,  and  VV.T.  Vetterling,  1980.  .Y amt  nail  Recipes.  Cam¬ 
bridge  University  Press. 

Robertson,  J.D.  and  D.  Corrigan,  1983.  R^uiation  patterns  of  a  shear-wave  vibrator  in  near-surface 
shale,  Geophysics ,  48,  19-26. 

Rokhlin,  S.  I.,  Bolland,  T.  K.,  and  Adler,  L.,  1986.  Reflection  and  refraction  of  elastic  waves  on  a 
plane  interface  between  two  generally  anisotropic  media,  J.  Acoust.  Soc.  Am,  70,  900  918. 

Saastamoinen,  P.R.,  1980.  On  propagators  and  scatterers  in  wave  problems  of  layered  elastic  media: 
a  spectral  approach.  Bull.  Seism.  Soc.  Am.,  70,  1125-1135. 

Schoenberg,  M.  and  Douma,  J.,  1988.  Elastic  wave  propagation  in  media  with  parallel  fractures 
and  aligned  cracks,  Geophys.  Prospecting.  36,  571-590. 

Shearer.  P.  M.  and  Chapman,  C.  IE,  1989.  Ray  tracing  in  azimuthally  anisotropic  media-I.  Results 
for  models  of  aligned  cracks  in  the  upper  crust,  Geophys.  J.,  96,  51-01. 

Smith,  13.  T.,  1970.  Matrix  Eigensystem  Routines  -EISPACK  Guide,  2nd  ed.,  Vol  0  of  Lecture 
Notes  in  Computer  Science  (New  York:  Springer- Velag). 

Stephen,  R.  A.,  19M.  Seismic  anisotropy  observed  in  upper  oceanic  crust.  Geophys.  IBs.  Lift..  S. 
*65-868. 

Takeuchi,  IE.  and  Saito.  M..  1972.  Seismic  surface  waves.  Method >  m  computational  Physic.-.  Aca 

ISO 


domic  Press,  New  York,  11,  217-295. 


Thomsen,  L.,  1986.  Weak  elastic  anisotropy,  Geophysics,  51,  1954-1966. 

Toksoz,  M.N.  and  H.H.  Kehrer,  1972.  Tectonic  strain  release  by  underground  nuclear  explosions 
and  its  effect  on  seismic  discrimination,  Geophys.  J.  R.  astr.  Soc.,  SI,  1  11  161. 

Wallace,  T.C.,  D.V.  Helmberger,  and  G.R.  Engen.  1983.  Evidence  of  tectonic  release  from  under¬ 
ground  nuclear  explosions  in  long-period  P  waves.  Bull.  Seism.  Soc.  Am..  IS.  593  613. 

Wilkinson,  J.  H.  and  Reinsch,  C.,  1971.  Linear  Algebra,  vol  II  of  Handbook  of  Automatic  Compu¬ 
tation,  Springer- Verlag,  New  York. 

Winterstein,  D.  F.,  1986.  Anisotropy  effects  in  P-wave  and  SH-wave  stacking  velocities  contain 
information  on  lithology,  Geophysics,  51,  661-672. 

Woodhouse,  J.  H.,  1974.  Surface  waves  in  a  laterally  varying  layered  structure,  Geophys.  J.  R.  astr. 
Soc,  37,  461-490. 


181 


APPENDIX  A. 

Wave  equation  with  out  body  force  (1) 

pd(U,i  —  T,jj 

—  L  Tij  —  L  Cijkl^kl  —  L  C ijklLufc 

( A-l) 

where 
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(A-5) 

Wave  equation  (e.g.  Kausel,  1986), 
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Now  stresses  in  equation  2  (main  text)  take  the  form 
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In  the  simpler  from 


-  ibJT  =  Vzua  4-  VyU,y  +  VZU'Z. 
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Substitute  the  trial  solutions,  we  will  get 


=  iu^Vr  +  PyVy)  uy  +VZ—  Uy 


( A- 15) 
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In  a  simplified  form 
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where 
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P  =  -V~x  and  0  =  PVZ 


Now  substitute  the  trial  solution  in  the  wave  equation  (A-8), 
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Note  that  Czz  =  Vz,  then  CZZP  —  -I  and  CzzO  =  -Vxy  The  /  is  3  x  3  identity  matrix.  After 
some  algebra  we  will  get  a  relation  between  TyZ  and  ( U,T ), 

Tz  =  iuj(QU  +  RT)  (A-20) 
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The  simplification  forms  of  0,Q  and  R  may  be  written  as 
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The  final  form  of  the  differential  system  will  be  expressed  as 
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which  is  the  equation  (3).  The  matrix  P  and  Q  are  symmetric. 


(A-25) 


187 


Figure  Captions 

Figure  1:  Anisotropic  layered  model. 

Figure  2:  Schematic  flowchart  of  the  algorithm. 

Figure  3:  Synthetic  seismograms  at  a  distance  of  75  km,  for  a  source  at  a  depth  of  10  km,  using 
the  results  of  Mandal  and  Mitchell  (1986)  (left  traces)  and  present  anisotropic  algorithm 
(right  traces).  Both  seismograms  are  computed  for  the  same  isotropic  model. 

Figure  4:  Schematic  diagram  of  an  anisotropic  half-space  model.  Dashed  lines  represent  frac¬ 
tures. 

Figure  5:  Three-component  synthetic  seismograms  at  a  5  km  epicentral  distance,  for  five  az¬ 
imuths  and  at  frequencies  of  1,  2,  3,  and  4  Hz.  The  anisotropic  model  has  5%  fracture  density. 
The  seismograms  at  4>  =  90°  are  equivalent  to  those  of  an  isotropic  reference  medium. 

Figure  0:  Three-component  synthetic  seismograms  at  a  50  km  epicentral  distance,  for  five  az¬ 
imuths  and  at  frequencies  of  1,  2,  3,  and  4  Hz  for  5%  fractured  model.  The  seismograms  at 
<t>  =  90°  are  equivalent  to  those  of  an  isotropic  reference  medium. 

Figure  7:  Three-component  synthetic  seismograms  at  a  5  km  epicentral  distance,  for  five  az¬ 
imuths  and  at  frequencies  of  1,  2,  3,  and  4  Hz.  The  anisotropic  medium  has  20%  fracture 
density.  The  seismograms  at  <f>  =  90°  are  equivalent  to  those  of  an  isotropic  reference  medium. 

Figure  8:  Three-component  synthetic  seismograms  at  a  50  km  epicentral  distance,  for  five  az¬ 
imuths  and  at  frequencies  of  1,  2,  3,  and  4  Hz  for  20%  fracture  density.  The  seismograms  at 
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4>  =  90°  are  equivalent  to  those  of  an  isotropic  reference  medium. 


Figure  9:  Schematic  diagram  of  an  anisotropic  layer  over  isotropic  half-space. 

Figure  10:  Three-component  synthetic  seismograms  for  a  layer  over  half-space  at  a  5  km  epicen- 
tral  distance,  for  five  azimuths  and  at  frequencies  of  1,  2,  3,  and  4  Hz.  The  seismograms  at 
4>  =  90°  are  equivalent  to  those  of  an  isotropic  reference  medium. 

Figure  11:  Three-component  synthetic  seismograms  for  a  layer  over  half-space  at  a  50  km  epi- 
central  distance,  for  five  azimuths  and  at  frequencies  of  1,  2,  3,  and  4  Hz.  The  seismograms 
at  4>  =  90°  are  equivalent  to  those  of  an  isotropic  reference  medium. 
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Figure  1:  Anisotropic  layered  model. 
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Kstep  Loop 

Layer  Loop 

Compute  Eigenvalues  and  Eigenvectors 
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Figure  2:  Schematic  flowchart  of  the  algorithm. 
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Isotropic  algorithm 


Anisotropic  algorithm 


— - 


Exp 1  o%  i on  source  depth  s  10  km 
Receiver  di*t*  =  75  km 


Figure  3:  Synthetic  seismograms  at  a  distance  of  75  km,  for  a  source  at  a  depth  of  10  km,  using 
the  results  of  Mandal  and  Mitchell  (1986)  (left  traces)  and  present  anisotropic  algorithm  (right 
traces).  Both  seismograms  are  computed  for  the  isotropic  model  given  in  figure  3. 


Figure  4:  Schematic  diagram  of  an  anisotropic  half-space  model.  Dashed  lines  represent  fractures. 
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Figure  5:  Three-component  synthetic  seismograms  at  a  5  km  epicentral  distance,  for  five  azimuths 
and  at  frequencies  of  1,  2,  3,  and  4  Hz.  The  anisotropic  model  has  5%  fracture  density.  The 
seismograms  at  <p  =  90°  are  equivalent  to  those  of  an  isotropic  reference  medium. 
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Figure  6:  Three-component  synthetic  seismograms  at  a  50  km  epicentral  distance,  for  five  azimuths 
and  at  frequencies  of  1,  2,  3,  and  4  Hz  for  5%  fractured  model.  The  seismograms  at  <p  =  90°  are 
equivalent  to  those  of  an  isotropic  reference  medium. 
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Figure  7:  Three-component  synthetic  seismograms  at  a  5  km  epicentral  distance,  for  five  azimuths 
and  at  frequencies  of  1,  2,  3,  and  4  Hz.  The  anisotropic  medium  has  20%  fracture  density.  The 
seismograms  at  <t>  =  90°  are  equivalent  to  those  of  an  isotropic  reference  medium. 
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Figure  8:  Three-component  synthetic  seismograms  at  a  50  km  epicentral  distance,  for  five  azimuths 
and  at  frequencies  of  1,  2,  3,  and  4  Hz  for  20%  fracture  density.  The  seismograms  at  <t>  =  90°  are 
equivalent  to  those  of  an  isotropic  reference  medium. 
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Figure  10:  Three-component  synthetic  seismograms  for  a  layer  over  half-space  at  a  5  km  epicen- 
tral  distance,  for  five  azimuths  and  at  frequencies  of  1,  2,  3,  and  4  Hz.  The  seismograms  at 
4>  -  90°  are  equivalent  to  those  of  an  isotropic  reference  medium. 
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Figure  11:  Three-component  synthetic  seismograms  for  a  layer  over  half-space  at  a  50  km  epicen- 
tral  distance,  for  five  azimuths  and  at  frequencies  of  1,  2,  3,  and  4  Hz.  The  seismograms  at 
<t>  =  90°  are  equivalent  to  those  of  an  isotropic  reference  medium. 


200 


FREQUENCY-DEPENDENT  ATTENUATION  IN  THE  CRUST 
M.  Nafi  Toksoz,  Batakrishna  Mandal  and  Anton  M.  Dainty 

Earth  Resources  Laboratory,  Massachusetts  Institute  of  Technology 

In  Preparation  for  Geophys.  Res.  Lett.,  October,  1989 


Abstract.  The  shear  wave  attenuation  ( Q -1)  in  the  crust  beneath  the  northeastern  US  and 
Scandinavia  is  determined  as  a  function  of  depth  and  frequency  with  the  aid  of  data  analysis  and 
theoretical  modeling  of  complete  wave  synthesis.  To  determine  attenuation  in  the  uppermost  crust, 
short  period  (0.25  <  T  <  2  sec)  Rg  waves  generated  by  explosive  sources  and  recorded  by  linear 
arrays  of  seismometers  are  used.  These  data  are  combined  with  longer  period  (6  <  T  <  20  sec) 
crustal  data  and  high  frequency  Lg  data  to  extend  the  investigation  into  the  deeper  crust.  The 
combination  of  broad-band  Rayleigh  waves  and  short  period  ( T  <  1.0  sec)  Lg  waves  makes  it 
possible  to  determine  attenuation  as  a  function  of  depth  and  frequency  in  the  crust,  parameterised 
as  Q  =  Qof  To  match  the  observations,  theoretical  models  for  three  mechanisms  that  contribute 
to  attenuation  in  the  crust — frequency  independent  anelastic  (Q~x)  and  frequency  dependent  fluid- 
flow  (QJ1)  and  scattering  (Qjl) — are  calculated.  The  best  fitting  models  require  low  Q  (Qo  <  100) 
in  the  shallow  crust  that  increases  with  increasing  frequency  as  f05  above  1  Hz.  In  the  middle  crust 
between  about  2  to  10  km  depth  Q  is  moderate  (100  <  Q o  <  500)  with  a  frequency  dependence 
of  f1.  In  the  lower  crust,  the  attenuation  is  very  low  ((J  >  500)  and  its  frequency  dependence  is 
difficult  to  establish.  The  low  intrinsic  Q  and  the  frequency  dependence  in  the  uppermost  crust 
can  be  explained  by  fluid  flow  in  fractures.  In  the  middle  crust  the  frequency  dependence  suggests 
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scattering.  The  increase  of  Q  in  the  lower  crust  is  most  likely  due  to  annealing  of  microcracks  under 
pressure  and  temperature  as  observed  in  the  laboratory. 

Introduction 

Seismic  attenuation  in  the  crust  and  lithosphere  affects  the  propagation  efficiency  of  seismic  phases 
at  regional  distances.  The  primary  objective  of  this  study  is  to  determine  shear  wave  attenuation 
(Q-1)  as  a  function  of  depth  and  frequency  in  the  crust.  The  anelastic  attenuation  of  seismic  waves  is 
commonly  determined  by  studying  the  decay  of  seismic  energy  with  epicentral  distances.  Different 
wave  types  (e.g.,  Lg,  Rg)  are  used.  Several  studies  have  evaluated  the  attenuation  characteristics  of 
seismic  waves,  especially  Lg  waves  in  various  parts  of  the  world  (e.g.  Nuttli,  1973;  Mitchell,  1975; 
Bollinger,  1979;  Nicolas  et  al. ,  1982;  Campillo  et  al. ,  1985).  This  subject  has  also  been  studied 
using  theoretical  approaches  (e.g.  Bath  and  Crampin,  1965;  Muller  and  Mueller,  1979;  Bache  et 
al.,  1981;  Bouchon,  1982;  Olsen  et  al.,  1983;  Herrmann  and  Kijko,  1983;  Campillo  et  al.,  1984). 
Mitchell  (1980,  1981)  has  evaluated  the  frequency  dependence  of  attenuation  of  Lg  waves  in  the 
United  States  for  longer  periods  (6  <  T  >  20  sec). 

In  this  article,  we  examine  the  attenuation  mechanisms  in  different  regions  of  the  earth’s  crust 
by  generating  complete  synthetic  seismograms  for  the  best  match  attenuation  model.  The  com¬ 
plete  synthetic  seismograms  are  computed  by  the  reflectivity  technique  combined  with  numerical 
wavenumber  integration  (Mandal  and  Mitchell,  1986;  Mandal,  1988).  The  effect  of  attenuation 
is  introduced  through  the  use  of  complex  velocities  (Aki  and  Richards,  1980,  Eq.  5.87,  p  182). 
We  assume  that  the  coda  and/or  Lg  attenuation  represents  the  intrinsic  attenuation  in  the  crust 
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(Toksoz  et  al.,  1988).  We  further  assume  that  three  mechanisms  contribute  to  intrinsic  attenuation: 
(1)  Frequency-independent  material  anelasticity;  (2)  frequency-dependent  fluid  flow  or  relaxation 
effects;  and  (3)  frequency-dependent  scattering.  Both  mechanisms  will  depend  on  depth  in  the 
crust.  For  layer  k ,  we  can  express  intrinsic  attenuation  as 

1  1  11 

Qtk(f)  ~  Q o*  +  Qfk(f)  Q,k(f)  (  } 

where  QtX(f)  is  the  total  attenuation,  Q"1  is  frequency-independent  anelastic  attenuation,  Qj1(f) 
is  the  frequency-dependent  attenuation  due  to  fluid  flow  or  relaxation,  and  is  the  frequency- 

dependant  attenuation  due  to  scattering.  Toksoz  and  Johnston  (1981)  summarise  results  for  anelas¬ 
tic  attenuation  in  dry  rocks,  demonstrating  that  it  is  frequency-independent.  Toksoz  et  al.  (1987) 
have  shown  that  P  and  S  wave  Q j  varies  as  fos  due  to  fluid  flow  above  a  critical  frequency  that 
may  be  as  low  as  1  Hz.  Dainty  (1981)  presents  a  model  of  scattering  attenuation  that  predicts  Qs 
varies  as  f1  and  suggests  that  it  should  apply  in  the  lithosphere  for  frequencies  greater  than  1  Hz. 
An  empirical  relationship  for  frequency  dependence  of  Q  is 

Qtk  =  Qok  /<•  (2) 

In  central  and  northeastern  U.S.,  eastern  Canada  and  European  shield  areas,  values  of  £  are  found 
to  vary  between  0.3  and  0.9  and  generally  cluster  around  £  =  0.5  (e.g.,  Singh  and  Herrmann, 
1983;  Pulli,  1984;  Campillo  et  al.,  1985;  Dysart  and  Pulli,  1987).  As  an  initial  attempt  to  fit  the 
data,  we  use  equation  (2)  and  then  interpret  the  value  of  £• 
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Q  Models 


Two  cases  of  modelling  for  different  observations  are  presented.  We  start  with  approximate  initial 
velocity  and  attenuation  models,  from  published  results  or  otherwise.  By  trial  and  error,  we 
establish  the  velocity  and  attenuation  models  of  the  area  which  match  the  observed  data  well. 
For  all  the  cases  we  use  the  short  period  vertical  component.  The  seismograms  are  observed  at 
ranges  5-250  km  and  are  in  the  frequency  band  0-10  Hz.  The  data  from  the  two  different  regions 
are: 

(a)  USGS  refraction  experiment  in  Maine  (Doll  et  al,  1986). 

(b)  Leningrad  quarry  blasts  observed  in  Fennoscandia  at  the  FINESA  array. 

USGS  refraction  experiment  in  Maine: 

The  seismograms  ( left  of  Figure  1)  from  6  to  29  km  were  acquired  in  northeast  Maine  (the  Chain 
Lakes  Massif)  during  the  1986  USGS  seismic  refraction  and  wide-angle  reflection  experiment  (Doll 
et  al .,  1986).  The  Rg  phase  is  well  recorded.  Figure  1  shows  the  comparison  between  the  observed 
seismograms  (left)  and  complete  wave  synthetics  (right).  To  model  this  region  we  start  with  an 
approximate  model  from  surface  wave  inversion  (Reiter  et  al,  1988).  The  final  velocity  and  atten¬ 
uation  model  are  shown  in  Figure  2.  The  source  is  assumed  to  be  a  point  vertical  force  at  depth 
40  m.  The  observed  and  synthetics  are  filtered  with  a  frequency  band  of  0.5  to  4  Hz  for  purposes 
of  comparison,  although  the  original  synthetics  were  calculated  for  0-10  Hz.  The  waveforms  are 
at  constant  scale  and  are  corrected  for  geometric  spreading.  Note  the  synthetics  do  not  generate 
scattered  waves  clearly  seen  in  the  observations.  This  suggests  that  the  real  earth  model  is  highly 
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heterogeneous  and  we  would  expect  some  scattering  attenuation.  In  the  Q  model,  we  see  £  =  0.5 
for  depths  less  than  1  km  increasing  to  (  =  1  below  this.  Below  8  km  depth  there  is  effectively  no 
constraint  on  Q. 

Leningrad  quarries: 

Two  quarry  blasts  near  Leningrad  observed  at  the  Fennoscandian  array  FINESA  are  presented 
here.  The  distances  are  200  and  250  km.  The  preliminary  velocity  model  was  taken  from  results 
from  the  FENNOLORA  experiment  published  by  Cassell  et  al.  (1983).  Observed  and  complete 
synthetic  seismograms  are  shown  in  Figure  3.  The  velocity  and  attenuation  models  used  for  the 
synthetic  seismograms  are  given  in  Figure  4.  The  quarry  blasts  are  represented  by  a  point  vertical 
force  at  depth  40  m  for  synthetic  computations.  Both  observed  and  synthetic  seismograms  are 
bandpass  filtered  between  0.8  and  4  Hz;  the  synthetics  were  originally  computed  for  a  0-8  Hz 
frequency  band.  Note  that  the  synthetics  are  not  corrected  for  instrument  response.  The  complete 
wave  synthetics  match  the  major  phases  (P,  Lg,  Rg)  of  the  observed  waveforms.  Specifically,  the 
synthetic  waveforms  predict  quite  accurately  the  higher  (Lg)  and  fundamental  (Rg)  modes.  Again, 
the  scattered  coda  is  not  predicted  by  the  synthetic  modelling.  There  is  a  large  phase  at  about  17 
sec  on  the  250  km  synthetic  which  is  e  reflection  off  the  Moho;  trial  calculations  demonstrate  that 
this  phase  is  very  sensitive  to  Moho  structure  and  may  be  easily  suppressed  by  slight  modifications 
of  it.  We  do  not  do  so  here  because  this  is  not  the  main  focus  of  the  study.  The  Q  model  again 
shows  a  thin  surface  layer  with  £  =  0.5  with  an  increase  to  (  =  1  between  2  and  16  km.  The  deep 
cruet  hae  a  high  shear  wave  Q  of  500  with  Q  =  0. 


DISCUSSION 


Both  of  these  cases  are  from  non-tectonic  regions  where  crystalline  rocks  are  found  at  the  surface. 
In  both  cases  the  seismic  attenuation  is  high  in  the  upper  crust  and  is  frequency  dependent.  Q0 
increases  with  depth  in  the  crust  from  less  than  100  in  the  top  2  km  to  100-500  from  2-10  km 
to  greater  than  500  in  the  lower  crust.  The  frequency  exponent  £  in  the  same  depth  ranges  is  0.5 
(0-2  km),  1.0  (2-10  km)  and  0  (lower  crust).  The  high  attenuation  and  the  frequency  exponent 
in  the  upper  2  km  can  be  explained  by  the  presence  of  fluid  flow  in  fractures.  Between  2  and  10 
km  the  frequency  exponent  indicates  that  scattering  is  the  dominant  mechanism  responsible  for 
attenuation.  In  the  lower  crust  anelastic  attenuation  is  occurring  with  the  high  Q  values  reflecting 
closing  of  fractures  and  annealing  of  microfractures  due  to  pressure  and  temperature  effects  (Smith 
and  Evans,  1984). 
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Figure  1:  Observed  and  synthetic  seismograms  as  a  function  of  distance;  observed  from  USGS 
refraction  experiment  in  Maine.  Synthetics  are  computed  for  a  point  vertical  force  of  depth 
40  m.  Receivers  are  on  the  surface  from  6-29  km  offset.  Both  observed  and  synthetics  are 
filtered  with  a  frequency  band  of  0.5-4  Hz. 

Figure  2:  Velocity  and  attenuation  models  for  synthetic  seismograms  in  Figure  1.  Solid  lines  refer 
to  P  wave  velocity  and  Q,  dashed  or  dotted  lines  to  S  waves.  For  Q,  Qa  -  2 Qg  and 

Q  =  Q0  x  for  both  P  and  S.  It  was  found  that  the  value  of  Qa  did  not  significantly 
affect  the  seismograms. 

Figure  3:  Comparison  of  observed  and  synthetic  seismograms  for  two  quarry  blasts  from  Leningrad 
recorded  at  the  FINESA  array.  Synthetics  and  observed  are  filtered  with  a  frequency  band 
of  0.8-4  Hz.  See  text  for  discussion. 

Figure  4:  Velocity  and  attenuation  models  for  synthetic  seismograms  in  Figure  3.  Symbols  as  for 
Figure  2. 
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Figure  2. 
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